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Abstract 


A  static-atomic  model  is  described,  which  may  be  employed  to  evaluate  the  characteristic 
vibrational  stiffness  of  an  atomic  lattice,  given  the  pair-wise  potential  of  the  constituent  atom. 
Because  the  vibrational  stiffness  is  directly  related  to  the  resultant  vibrational  frequency 
spectrum  of  the  lattice,  the  method  may  be  used  to  infer  the  behavior  of  the  characteristic  lattice 
frequency  as  a  function  of  lattice  spacing.  The  characteristic  frequency  behavior  is  sufficient  to 
determine  the  Griineisen  function,  an  important  thermodynamic  parameter  relating  to  thermal 
behavior  of  a  crystal  lattice.  The  current  method  computes  and  utilizes  several  spring  constants 
derived  from  a  static  lattice  in  order  to  infer  the  characteristic  vibrational  behavior.  No  atomic 
dynamics  calculations  involving  either  the  equations  of  motion  or  modal  (vibrational)  analysis 
are  required.  As  such,  the  method  generally  requires  mere  seconds  of  computation  on  today's 
generation  of  desktop  workstations.  Results  indicate  that  the  vibrational  stiffness  of  the  lattice 
is  qualitatively  distinct  from  the  volumetric  stiffness  of  the  lattice,  and,  furthermore,  that  the 
resulting  lattice  behavior  can  be  described,  over  a  wide  region  of  lattice  spacing,  by  an  analytical 
equation  of  state  in  terms  of  lattice  frequency. 
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PART  I — A  Static  Atomic  Paradigm  for  Lattice  Vibration 


1.  Introduction 

The  determination  of  the  vibrational  behavior  of  a  crystal  lattice  is  fundamental  to  the 
understanding  of  the  thermal  behavior  of  crystalline  materials.  The  magnitude  and  distribution 
of  the  lattice’s  vibrational  spectrum  was,  for  example,  a  key  ingredient  in  the  development  of 
plausible  specific-heat  models  by  Debye;  Bom  and  von  Karman  [1];  and  others  early  this  century. 
Whereas  the  specific  heat,  Cv,  relates  the  behavior  of  thermal  energy  and  temperature  at  constant 
volume,  V,  it  is  the  Griineisen  function,  T,  macroscopically  given  by  T  =  -V  (dpth  ldEth)v,  that 
relates  the  behavior  of  a  lattice’s  thermal  energy,  Elh,  to  the  thermal  pressure,  pth,  at  constant 
volume. 

At  the  atomic  level,  there  is  contribution  to  the  Griineisen  function,  associated  with 
each  mode  of  the  lattice’s  vibrational  spectrum,  which  is  given  by  r,  =  -V(dv,  /3V)T/v,  ,  where 
V,  denotes  a  frequency  component  of  the  vibrational  spectrum.  Historically  [2],  attempts  to 
aggregate  the  T,  components  into  a  macroscopic  value  for  T  were  accomplished  with  the  use  of 
two  assumptions:  (1)  that  the  frequency  spectrum  took  on  a  particular  idealized  form  (e.g.,  that 
of  Debye);  and  (2)  that  the  behavior  of  every  v,  component  changed  with  volume  in  a 
characteristically  similar  manner.  With  these  assumptions,  knowledge  of  the  behavior  of  the 
lattice’s  characteristic  frequency  alone  (and  not  the  details  of  the  vibrational  spectrum)  is 
sufficient  to  determine  the  macroscopic  value  for  T,  as  T  =  -V  (d(o/5V)r  /to,  where  to  is  the 
characteristic  frequency  of  the  lattice,  rather  than  a  constituent  component  of  the  spectrum.  With 
the  use  of  one  further  assumption,  generally  supported  by  data,  that  T  is  independent  of 
temperature,  T,  the  definition  for  the  Griineisen  function  may  be  alternately  given  as 

l/\|f  =  T/V  =  -(dto/JV)/co  ,  (1-1) 

where  \|t  is  introduced  in  this  equation  as  a  definition.  Since  0)  is  constant  (dto/dV  is  zero)  for 
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the  special  case  of  a  harmonic  lattice,  T  may  also  be  thought  of  as  a  measure  of  the  lattice’s 
vibrational  anharmonicity. 

More  recently,  Plendl  [3, 4]  analyzed  actual,  rather  than  idealized,  frequency  spectra. 
Through  extensive  correlation,  he  quantitatively  determined  that  the  characteristic  frequency  of 
any  given  atomic  vibrational  spectrum  could  be  computed  by  the  center  of  gravity  of  that 
spectrum,  the  so-called  centro-frequency,  as  in 


jv/(v)  dx 

«  =  V  -  ± -  .  d-2) 

j/(v)  dx 


where  xctr  is  the  centro-frequency,  and  /(v)  is  the  experimentally  measured  frequency  spectrum. 
That  the  characteristic  frequency,  co,  may  be  given  by  xctr  is  telling.  It  says  that  the  actual 
spectrum  must  straddle  the  characteristic  frequency,  unlike,  for  example,  the  idealized  Debye 
spectrum  in  which  the  characteristic  frequency  is  the  largest  frequency  component  of  the 
spectrum.  More  specifically,  the  characteristic  frequency  is  a  weighted  average  of  the  frequency 
spectrum.  This  concept  is  of  some  importance  to  subsequent  discussions  on  the  nature  of  the 
characteristic  lattice  stiffness. 


2.  Previous  Models  of  the  Griineisen  Function 

A  number  of  models  have  been  introduced  for  the  calculation  of  T.  An  early  model,  of 
Slater  [2],  examined  the  restoring  force  upon  an  atom  displaced  from  its  rest  position  in  the  lattice 
and  used  this  to  determine  a  spring  constant  and  thus  a  frequency.  However,  because  Slater 
quantified  the  spring  constant  in  terms  of  pressure-volume  compressibility,  rather  than  an 
interatomic  force-displacement  relation,  his  formula  could  not  recreate  the  proper  value  of  T  =  0 
for  the  limiting  harmonic  case.  Dugdale  and  MacDonald  [5]  remedied  this  deficiency  by  using 
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the  properly  constructed  spring  constant.  However,  in  both  of  these  analyses,  only  nearest- 
neighbor  interactions  were  considered  and  only  longitudinal  waves  were  modeled. 

Later,  Pastine  [6]  and  Vashchenko  and  Zubarev  [7]  introduced  separate  models  for  T, 
which  tried  to  improve  on  prior  work  by  accounting  for  transverse  waves,  in  addition  to  the 
longitudinal  waves  considered  in  earlier  works.  Both  these  works  relied  on  the  characterization 
of  empirically  composed,  volume-dependent  Poisson  ratios,  the  theoretical  basis  for  which  was 
no  better  understood  than  that  for  I\  unfortunately.  In  the  case  of  Vashchenko  and  Zubarev, 
however,  the  empirical  function  was  chosen  in  such  a  manner  as  to  reproduce,  for  three  special 
cases,  the  models  of  Slater,  Dugdale  and  MacDonald,  and  the  free  volume  theory,  respectively. 

Common  to  these  and  other  analyses  [2, 5, 8, 9]  was  the  assumption  that  the  volumetric 
stiffness  of  the  lattice,  proportional  to  the  second  derivative  of  the  energy  potential  with  respect 
to  lattice  spacing,  d2Epot/dX2,  is  proportional  to  the  vibrational  stiffness  of  the  lattice.  Such  an 
assumption  has  been  pervasive  for  two  reasons:  (1)  for  a  harmonic  lattice,  it  is  strictly  true;  and 
(2)  when  only  nearest-neighbor  interactions  are  considered  for  the  case  of  infinitesimal  vibrational 
amplitudes,  it  is  true  by  virtue  of  linearization.  These  two  assumptions  cover  a  wide  range  of 
modeling  strategies,  as  the  only  alternative  would  seem  to  be  the  modeling  of  nonlinear  effects 
and  non-nearest-neighbor  interactions,  with  the  resultant  set  of  coupled,  nonlinear  equations  being 
solvable  through  numerical  means  only  ( e.g .,  time  integration  or  modal  analysis). 

Recently,  Segletes  [10]  proposed  a  semi-empirical  equation  of  state  that  accurately 
describes  the  lattice  force  and  vibrational  stiffness  for  many  metals  to  several  megabars  of 
pressure,  in  terms  of  the  vibrational  frequency  of  the  lattice.  In  this  model,  the  zero-temperature 
lattice  force  varies  as 


F  =  (C0/ro)2  -Z2/v  •  (co/co0f  ln(w/co0)  ,  (2-1) 

where  the  “0”  subscript  denotes  conditions  at  the  reference,  zero-force  state,  C0  is  the  bulk 
sound  speed  at  the  reference  state,  and  K,  taken  as  a  constant  K  =  r|/(3r0),  indicates  a  ratio  of 
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volumetric  to  vibrational  anharmonicity  (where  tj  is  a  nondimensional  parameter  employed  by 
Rose  et  al.  [11]),  which  correlates  the  curvature  of  the  lattice  potential  (i.e.,  the  slope  of  the 
force)  at  the  reference  state  to  the  binding  energy  of  the  lattice  potential.  The  variable  k  is  the 
current  lattice  spacing,  and  \|r,  as  defined  in  eqn  (1-1),  is  simply  the  ratio  V/T.  The  lattice  force, 
F,  is  characterized  by  pc  k2,  where  pc  is  the  zero-temperature  (so-called  cold)  pressure  in  the 
lattice.  Subsequent  analysis  of  this  equation-of-state  form  [12-14]  revealed  that  the  volumetric 
stiffness  of  the  lattice,  denoted  dF/dk,  was  of  a  different  functional  form  than  the  vibrational 
stiffness  of  the  lattice,  denoted  dFvlb  /dx.  In  particular,  for  vibrations  of  infinitesimal  amplitude 
(which  we  assume  for  the  low  temperature  behavior  of  solids),  the  vibrational  stiffness,  from 
considerations  of  a  linearized  vibration  analysis,  varies  as 

-dFvib/d x~(D2  .  (2-2) 

But  rather  than  the  volumetric  stiffness  following  this  same  functional  form,  a  more  complicated 
expression  arises  from  eqn  (2-1),  which,  even  when  idealized  to  a  quasi-harmonic  form,  varies 
instead  along  the  lines  of 


-dF/dk  -  to  [  ln(to/w0)  +  1  ]  (quasi-harmonic).  (2-3) 

This  disparity  provoked  speculation  that:  (1)  the  omission  of  non-nearest-neighbor  interactions 
from  the  early  analyses  may  have  led  to  a  systematic  error,  and  (2)  there  might  be  some  hope 
of  capturing  the  interactions  of  nonnearest  neighbors  within  the  context  of  an  analytical  form 
along  the  lines  of  eqn  (2-1). 

It  is  the  primary  intent  of  this  report  to  quantitatively  analyze  the  influence  of  non-nearest- 
neighbor  interactions  on  various  stiffnesses  of  the  lattice,  in  hopes  of  ascertaining  whether  or  not 
there  is  any  basis  for  the  distinction  cited  between  eqns.  (2-2)  and  (2-3),  which  might  lead  to  a 
relation  like  eqn  (2-1).  To  compute  the  details  of  the  vibrational  spectrum,  a  modal  analysis  or 
time-integrated  dynamic  solution  would  need  to  be  performed  on  a  large  array  of  lattice  atoms. 
But  since  only  the  characteristic  frequency  is  desired,  which  is  sufficent  to  obtain  a  value  for  T, 
we  seek  a  shortcut  to  the  full  dynamic  solution.  Unlike  Slater  [2]  and  Dugdale  and 
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MacDonald  [5],  the  current  analysis  considers  the  behavior  of  both  the  longitudinal  and  transverse 
stiffnesses.  However,  it  does  so  along  only  the  primary  axes  of  the  crystal.  A  contemplated 
future  effort  will  work  to  address  this  limitation.  Finally,  the  method  to  be  proposed  here  uses 
a  given  interatomic  pairwise  potential  (energy  vs.  distance)  in  order  to  infer  the  characteristic 
frequency  behavior  of  the  lattice.  It  therefore  goes  without  saying  that  the  validity  of  any  results 
obtained  here  is  bounded  by  the  domain  of  applicability  of  the  selected  potential.  We  intend  to 
show  that,  in  spite  of  ignoring  directional  variations  in  crystalline  behavior,  a  meaningful 
description  of  the  characteristic  frequency  of  the  lattice  can,  nonetheless,  be  obtained. 

3.  Non-Nearest-Neighbor  Interactions 

One  may  argue,  for  several  reasons,  that  the  vibrational  and  volumetric  stiffnesses  of  an 
atomic  lattice  are  not  functionally  equivalent.  One  compelling  argument  [13, 14]  may  be  made 
by  considering  a  one-dimensional  string  of  atoms.  The  volumetric  stiffness  quantifies  the 
resultant  change  of  externally  applied  force  required  to  alter  the  spacing  between  every  adjacent 
atom  in  the  lattice  by  an  amount  dk,  whereas  the  vibrational  stiffness  quantifies  the  restoring 
force  on  a  single  atom  in  question  when  it,  as  part  of  a  moving  wave,  is  displaced  a  distance  dx 
with  respect  to  a  lattice  that  otherwise  remains,  on  average,  at  a  fixed  location.  In  the  former 
case,  the  resultant  change  in  separation  to  the  n’th  nearest  neighbor  is  n  dk,  while  in  the  latter 
case,  it  is  merely  dx.  Clearly,  these  two  situations,  and  thus  the  two  associated  stiffnesses,  must 
be  different  if  interactions  are  considered  for  neighbors  n>  1.  Furthermore,  Plendl[3, 4]  has 
actually  reported  experimental  data  on  this  distinction,  by  comparing  his  centro-frequency  (a 
measure  of  vibrational  stiffness)  to  the  so-called  “definite  frequency”  (a  measure  of  volumetric 
stiffness)  at  the  equilibrium  states  for  many  ionic  crystals.  He  noted  that  the  interrelation  of 
these  two  frequencies  could  be  characterized  in  terms  of  the  lattice’s  exponent  of  repulsion  and 
they  were  observed  to  vary  from  each  other  by  as  much  as  ±20%. 

If  one  considers  a  single,  infinitesimal  wave  of  disturbance  along  the  principal  axis  of  a 
crystal  lattice,  then  the  characteristic  vibrational  stiffness  of  the  lattice  may  be  acquired  by 
finding  the  restoring  force  on  one  atom  of  interest  (atom  O)  by  moving  a  plane  of  atoms 
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containing  O  a  distance  dx  with  respect  to  an  otherwise  stationary  lattice.  If  the  displacement, 
dx,  is  perpendicular  to  the  plane  of  atoms,  the  vibration  is  a  longitudinal  one.  If  dx  is  parallel 
with  the  plane,  the  vibration  is  transverse.  If  knowledge  of  the  characteristic  stiffness  behavior 
can  be  obtained  in  this  manner,  then  vibration  theory  for  infinitesimal  vibrations  [i.e.,  eqn  (2-2)] 
will  immediately  yield  the  characteristic  vibrational  frequency  behavior. 

For  the  more  realistic  case  of  finite  vibrational  amplitudes  and  multiwave  disturbances, 
the  instantaneous  interatomic  stiffness  will  vary  finitely  in  time  because  the  instantaneous 
separation  of  atoms  will  be  a  function  of  the  frequency  phase-shift  of  adjacent,  vibrating  atoms. 
The  averaged  value  of  this  time-dependent  stiffness  may,  thus,  no  longer  be  exactly  that  of  the 
infinitesimal-amplitude  single-wave  stiffness  described  in  the  proceeding  text,  though  we  hope 
the  infinitesimal  stiffness  will  still  serve  as  a  reasonable  approximation.  In  a  dynamic  solution, 
the  instantaneous  distance  between  any  two  atoms  (i.e.,  phase  information),  because  of  vibration, 
will  be  sometimes  larger  than  and  other  times  smaller  than  the  current  equilibrium  rest  distance. 
The  instantaneous,  pairwise  stiffness  will  also  fluctuate  about  this  equilibrium  value  (unless  the 
interaction  is  harmonic,  in  which  case  the  stiffness  is  independent  of  separation).  During  the 
course  of  the  vibration,  as  the  stiffness  is  locally  larger  than  the  current  equilibrium  stiffness, 
higher  frequency  modes  tend  to  be  excited,  and  correspondingly,  when  the  instantaneous  stiffness 
is  lower  than  equilibrium,  lower  frequency  modes  are  excited.  Since  a  vibration  necessarily 
makes  excursions  about  both  sides  of  the  equilibrium  point,  both  higher  and  lower  frequency 
modes  will  be  excited  when  anharmonicity  exists,  as  compared  with  a  locally  harmonic  spectrum. 

And  yet,  the  average  effect  on  the  centro-frequency  of  the  spectrum  may  be  small  by 
comparison,  since  the  simultaneous  introduction  of  higher  and  lower  frequency  modes  will  tend 
to  cancel  out  in  the  centro-frequency  calculation  of  eqn  (1-2).  So  it  goes,  too,  with  the 
vibrational  stiffness  because  of  the  monotonic  relation  between  frequency  and  stiffness,  eqn  (2-2). 
Thus,  while  anharmonic  effects  will  cause  the  pairwise  stiffness  to  fluctuate  in  time  about  its 
equilibrium  value  for  finitely  large  vibrational  amplitudes,  it  is  the  static  value  of  that  stiffness 
that  might  reasonably  be  used  to  approximate  the  time-averaged  characteristic  stiffness 
corresponding  to  the  centro-frequency. 
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In  the  same  way  that  determination  of  a  spring  stiffness,  for  a  simple  mass-spring 
arrangement,  is  sufficient  to  characterize  the  resultant  vibrational  frequency  without  having  to 
reevaluate  the  equations  of  dynamic  motion,  we  intend  to  compute  and  utilize  several  spring 
constants  derived  from  a  static  lattice  in  order  to  infer  the  characteristic  vibrational  behavior.  No 
atomic  dynamics  calculations  involving  either  the  equations  of  motion  or  modal  analysis  will  be 
required.  Essentially,  the  equation  of  motion  is  already  contained  in  the  solution  represented  by 
eqn  (2-2),  for  a  simple  spring-mass  system  under  infinitesimal  vibrational  amplitude. 

One  other  reason  to  believe  that  the  static  stiffness  of  a  lattice  may  be  used  to  characterize 
the  dynamic  solution  is  that,  in  a  lattice,  vibrational  amplitude  is  a  measure  of  temperature.  And 
since  T  is,  for  actual  lattices,  generally  found  to  be  temperature-insensitive  over  a  large  range, 
we  may  conclude  that  the  stiffness  and  its  derivatives  are  vibrational-amplitude  insensitive  to  the 
same  extent. 

4.  Volumetric  Distention 

Let  8(5)  denote  the  energy  potential  between  any  two  atoms  in  a  3-D  simple-cubic  lattice, 
separated  by  a  distance  s.  Then  /,  given  by  -de/ds  (positive  in  compression),  will  denote  the 
pairwise  force.  Considering  our  atom  of  interest,  atom  O,  to  be  at  the  origin,  we  use  the  index 
ijk  to  denote  an  atom  located  relative  to  the  origin  at  coordinate  (x,  y,  z)  =  ( Xi,Xj,Xk ),  where  X 
is  the  current  lattice  spacing  of  the  given  simple-cubic  lattice.  Applying  the  distance  formula, 
sijk  =  C*2  +  y2  +Z2),  from  our  atom  of  interest  to  the  atom  ijk  is  simply 

S'^xie+f+k2)'12  .  (4-1) 

If  we  use  Epol(k )  to  denote  the  summed  potential  energy  for  the  atom  at  the  origin,  arising  from 
the  contributions  of  every  other  atom  in  the  lattice,  then  for  a  lattice  of  infinite  extent  with  an 
instantaneous  lattice-spacing  of  X,  we  have 
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(4-2) 


oo  oo  oo 

■  E  E  E  <=<V 

l--oo  y-  -oo  k  -  -oc 


where  we  define  e(0)  =  0,  to  preclude  the  situation  of  atom  O  producing  a  contribution  to  its  own 
potential. 

If  we  first  consider  the  case  of  volumetric  expansion,  whereby  the  spacing  between 
adjacent  atoms  in  the  lattice  is  increased  from  X  to  X  +  dX,  eqn  (4-1)  indicates,  for  a  given  pair 
of  atoms,  that  sijk  is  directly  proportional  to  X.  Thus,  the  separation  metric  between  the  origin 
and  atom  ijk,  for  volumetric  expansion,  is  obtained  from  eqn  (4-1)  as 

dsvol/dX  =  sljk/X  =  (i2+f  +  k2)v2  .  (4-3) 

Eqn  (4-3)  merely  indicates  that  the  amount  of  increased  separation  caused  by  an  incremental 
change  in  X  is  directly  proportional  to  the  distance  of  ijk  from  O.  The  quantity  dsml/dX  is  a 
function  only  of  the  topological  indices,  i,  j,  and  k,  and  not  of  X  or  s.  Thus,  for  any  given  atom, 
ijk,  this  metric  is  constant  under  volumetric  dilatation.  As  such,  higher  derivatives,  such  as 
d2svol/dX2,  are  identically  zero. 

A  spring  force,  in  general,  is  given  by  F  ~  dE/dx.  By  attaching  to  the  term  “lattice  force” 
the  physical  meaning  of  pc  X2,  where  pc  is  the  cold  pressure  of  the  lattice,  it  may  be  shown,  for 
the  current  case,  that  F  =  ~(dEpol/dX)/3,  consistent  with  the  spring  notion.  We  therefore  compute 
F  (physically  representing  the  external  unit  force  required  to  volumetrically  compress  the  lattice) 
as 


F(X)  =  -1/3  £  £  £  deJdX\s  .  (4-4) 

I  -  -oo  j-~oo  fc-  ~ 
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Furthermore,  the  term  dddX,  for  the  case  of  volumetric  expansion,  is  ddds  -dsvol/dX.  The  first 
deriviative  of  the  pairwise  energy  potential  is  simply  the  negative  of  the  pairwise  force — thus, 


m  =  1/3  e  E  E  /(»„,)  • 

j--oo  -oo  lc-~°o 


(4-5) 


We  wish  to  compute  the  spacial  derivative  of  eqn  (4-5)  for  the  case  of  volumetric  deformation, 
since  the  force-derivative  is  itself  one  of  the  spring  stiffnesses  being  sought.  For  the  case  of 
volumetric  distention,  we  again  take  the  derivative  with  respect  to  A.  to  obtain 


dF/dk  =  1/3  E  E  E  df/ds  L,  WJ^lvk)2  •  (4-6) 

/‘sr-OO  j—  -  OO  k-  "OO 


With  the  use  of  eqn  (4-3),  eqns  (4-5)  and  (4-6)  may  be  simplified  as 


F  =  1/3  E  E  E  /(^y*)  (i2+/  +  ^2)1/2  (4-7) 

/=-oo  j=~cc  k~-oo 


and 


=  1/3  E  E  E 

/'=-oo  j~ -oo  k=-°o 


(i2  +j2  +  k2) 


(4-8) 


5.  Vibrational  Distortion  (Longitudinal) 


In  contrast  to  the  situation  of  volumetric  distention,  an  instantaneous  snapshot  of  the 
simplest  vibrational  distortion  ( i.e a  planar  pulse)  involves  the  displacement  of  a  single  plane 
of  atoms  with  respect  to  an  otherwise  stationary  lattice.  Again,  consider  our  atom  of  interest,  O, 
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to  be  initially  located  at  the  origin,  at  the  moment  a  longitudinal  plane  wave  traveling  from  the 
+x  direction  to  the  —x  direction  just  passes  the  origin.  At  this  instant,  all  atoms  in  the  x  =  0  plane 
are  displaced  a  distance  dx  in  the  negative  x  direction,  whereas  every  other  atom  in  the  lattice 
is  located  at  its  nominal,  undisturbed  position.  For  all  atoms  in  the  displaced  plane,  no  change 
in  position  occurs  with  respect  to  each  other,  since  all  atoms  in  that  plane  are  moving  in  unison. 
However,  we  may  compute  the  vibrational  separation  metric  between  atom  O  and  atom  ijk  not 
located  in  the  x  =  0  plane.  In  terms  of  the  motion  of  atom  O  in  the  x  direction,  we  can  quantify 
the  vibrational  separation  metric  by  taking  the  derivative  of  the  atomic  separation  to  atom  ijk, 
sifk  =  C*2  +  y2  +Z2),  with  respect  to  x  only.  We  obtain 

dsvib  /dx  =  x/sijk  =  cos  y  (/  *  0  for  longitudinal  vibrations),  (5-1) 

where  cosy  is  the  direction  cosine  between  the  line  of  action,  from  O  to  ijk,  and  the  vibrational 
movement  along  the  x  axis.  For  this  particular  vibrational  mode,  it  is  given  by 
cos  y  =  z‘/(z2  +/  +  k2)112.  The  i  *  0  restriction  in  eqn  (5-1)  accounts  for  the  fact  that  all  other  atoms 
in  the  x  =  0  plane  are  moving  in  unison  with  atom  O,  and  their  vibrational  separation  metric  thus 
remains  zero.  Unlike  dsvol  IdX,  the  vibrational  separation  metric  is  not  independent  of  changes 
in  x.  Therefore,  higher  derivatives  will  need  to  be  computed  as  well;  for  example, 

d2sub  /dx2  =  1  lsijk  -  x1/s,jik  =  (sin2  y)ls,jk  (z  *  0),  (5-2) 

and  so  on.  As  a  point  of  note,  for  the  vibrational  mode  being  considered, 
sin  y  =  [(/+  k2)/(z'2  +/  +  k2)]1'2. 


The  longitudinal  vibrational  force  in  the  lattice,  call  it  FL  (physically  representing  the 
resultant  force  on  the  oscillating  atom,  O),  and  its  derivative  may  be  obtained  in  a  like  manner 
to  the  volumetric  force  quantities,  though  they  will  differ  from  the  earlier  results  in  that  the 
derivative  of  e  is  taken  with  respect  to  the  vibrational  displacement,  x,  and  not  the  lattice  spacing, 
X.  The  vibrational  force  and  its  derivative  are  thus 
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FlO^>  =  E  E  E  dSJdX\ijk  (i  *  °)  ’  (5“3) 


and 


dFt/dx  =  Y,  Y  E  [  1^  (d^/dxi,..*)2  +  /(^)  avfe/a*2  ]  o'*o)  .  (5-4) 


Because  the  sum  is  performed  at  the  moment  when  the  vibrating  atom  is  instantaneously  at  its 
equilibrium  position  (the  origin),  the  resultant  force,  Fvih,  on  the  atom  should  be  momentarily 
zero.  We  see  from  eqn  (5-1)  that  Bsvib  /dx  in  the  -x  half-space  is  opposite  in  value  to  the 
corresponding  point  in  the  +x  half  space.  Thus,  the  sum  of  eqn  (5-3),  when  taken  over  the  full 
3-D  space,  will  indeed  be  zero.  Were  the  sum  to  be  performed  over  the  +x  half-space  only,  the 
resulting  nonzero  value  for  force  would  represent  the  local  force  being  exerted  on  the  +x  side  of 
the  atom.  However,  though  the  net  force  on  an  unextended  spring  is  zero,  the  stiffness  of  that 
same  spring  is  not  zero.  Eqn  (5-4)  reveals  this  point  to  be  the  case,  as  the  squared  term  will 
always  yield  a  positive  contribution  to  the  summation.  It  is  with  the  use  of  eqns  (5-1)  and  (5-2) 
that  eqns  (5-3)  and  (5-4)  may  be  simplified  as 


fl  -  £  £  £  /<  v  cosv  "  0  (>■  *  0)  .  <5-5> 


and 


dFJdx  =  Y  E  E  [  dfids  \sj  cos2y  +  f(stjk)/s.jk-  sin2y  ]  (i  *  0)  .  (5-6) 

f=-oe  j~~oo  k~ 
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Eqn  (5-6)  confirms  what  we  know  to  be  true,  were  the  lattice  to  be  composed  of  harmonic  (i.e., 
linear)  springs.  That  is,  when  the  force-displacement  relationship  for  the  constituent  springs  of 
a  lattice  is  linear,  such  that  dfds  =  fls,  the  characteristic  vibrational  stiffness  of  the  lattice  is 
proportional  to  the  stiffness  of  the  component  spring.  Furthermore,  that  constant  of 
proportionality  is  the  number  of  springs  acting  on  the  atom. 

6.  Vibrational  Distortion  (Transverse) 

Like  the  longitudinal  vibration,  an  instantaneous  snapshot  of  a  transverse  vibrational 
distortion  also  involves  the  displacement  of  a  single  plane  of  atoms  with  respect  to  an  otherwise 
stationary  lattice.  Again,  we  consider  our  atom  of  interest,  O,  to  be  initially  located  at  the  origin, 
at  the  moment  a  transverse  plane  wave  traveling  from  the  +>•  direction  to  the  —y  direction  just 
passes  the  origin.  At  this  instant,  all  atoms  in  the  y  =  0  plane  are  displaced  a  distance  dx  in  the 
negative  x  direction,  whereas  every  other  atom  in  the  lattice  is  located  at  its  nominal,  undisturbed 
position.  For  all  atoms  in  the  displaced  plane,  no  change  in  position  occurs  with  respect  to  each 
other,  since  all  atoms  in  that  plane  are  moving  in  unison.  However,  we  may  compute  the 
vibrational  separation  metric  between  atom  O  and  atom  ijk  not  located  in  the  y  =  0  plane.  The 
resulting  separation  metric  between  O  and  atom  ijk,  when  our  atom  of  interest  is  momentarily 
displaced  to  coordinate  (-d!x,0,0),  is  likewise  obtained,  in  the  manner  of  eqn  (5-1),  as 

dsvib/dx  =  x/sljk  =  cosy  (y  *0  for  transverse  vibrations)  .  (6-1) 

The  transverse-vibrational  force  in  the  lattice,  call  it  FT  (physically  representing  the 
resultant  force  on  the  oscillating  atom,  O),  and  its  derivative  may  be  obtained  in  a  like  manner 
to  the  longitudinal-vibrational  force  quantities.  The  transverse-vibrational  force  and  its  derivative 
are  thus 


FA)  -  E  E 

I  -  -oo  y-  -oc 


E  dsJdx  * 


0*0)  , 


(6-2) 
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and 


dFT/dx  =  £  E  £  t  l,w  (d*v,*/a*  ijk )2  +  /(^)  ]  O'  *  °)  •  (6~3) 

/  -oo  j~-rx>  A'=-oo 


Like  the  longitudinal  case,  the  resultant  force,  FT,  on  the  atom  will  be  momentarily  zero,  though 
the  stiffness  of  that  same  spring  is  not  zero.  We  see  that  the  only  difference  between  these  and 
the  longitudinal  equations  is  the  nature  of  the  summation  restriction  as  to  which  atoms  are 
moving  in  unison  with  atom  O.  The  final  results  are  thus 


oo  oo  oo 

Ft  =  £  £  £  /Oy*)  COSY  =0  0  *0),  (6-4) 

i~-  oo  j=- oo  k-~ 


and 


dFT/dx  =  £  £  £  [  1^  -  cos2y  +  f(stjk)/sijk-  sin2Y  ]  0*0)  •  (6-5) 

i~~oo  oo  k--oo 


One  last  note  regarding  transverse  vibrations  is  that  transverse  waves  traveling  in  the  y  and  z 
directions  will  both  produce  displacements  in  the  x  direction,  whereas  only  a  longitudinal 
vibration  in  the  x  direction  will  do  so.  Thus,  there  are  two  possible  transverse  vibrational 
components  for  every  longitudinal  component.  Because  of  considerations  of  symmetry,  however, 
both  transverse  wave  components  will  exhibit  characteristic  behavior  identical  to  that  expressed 
by  eqn  (6-5). 

7.  Vibrational  Stiffness  (Aggregated)  and  Its  Derivatives 

The  vibrational  spectrum  of  an  actual  material  is  composed  of  many  different  modes.  It 
is  plain  to  see  that  equations  developed  to  this  point,  for  longitudinal  and  transverse  vibrations, 
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consider  only  the  modes  along  the  principal  axes  of  the  lattice.  In  this  current  work,  we  limit 
the  analysis  to  this  restriction,  and  for  two  reasons:  (1)  primarily,  for  the  sake  of  simplicity;  and 
(2)  because  past  works  by  Slater  [2]  and  Griineisen  suggest  that  different  vibrational  modes  vary 
in  a  characteristically  similar  manner  under  changes  in  lattice  spacing.  Though  not  reported  here, 
we  have  formulated  modified  relations  to  examine  other  vibrational  modes  of  the  lattice. 
Preliminary  results  do  indicate,  to  a  large  extent,  a  characteristically  similar  behavior  in  terms 
of  variation  of  the  modal  frequency  with  lattice  spacing.  Ideally,  we  hope  in  a  future  effort  to 
integrate  analytically  all  of  the  vibrational  modes  into  a  characteristic  stiffness. 

Even  with  the  simplification  of  restricting  the  study  to  the  principal  vibrational  modes  of 
the  lattice,  we  still  have,  at  this  point,  a  need  to  aggregate  the  one  longitudinal  and  two 
corresponding  transverse  vibrational  stiffnesses  into  a  net  vibrational  stiffness  quantity.  It  may 
be  shown  from  the  work  of  Slater  [2]  or  Brillouin  [8]  (which  presents  the  result  directly)  that  the 
frequency  components  are  related  to  an  aggregated  frequency  by  way  of 

3/g)3  =  1  /fa)3  +  2/ wl  .  (7-1) 

These  frequencies  vary  with  their  associated  stiffnesses  by  way  of  eqn  (2-2).  Therefore,  an 
aggregated  stiffness  may  be  obtained  as 

-dFJdx  =  <3  /  [1  /{-dFJdxf2  -  2/(-dFT/dx)y2]  .  (7-2) 

Note  that  even  though  the  longitudinal  and  transverse  components  of  the  vibrational  stiffness  may 
change  sign  (becoming  positive)  at  large  extensions  of  the  lattice,  eqn  (7-2)  loses  its  traditional 
meaning  when  either  of  the  component  stiffnesses  becomes  positive,  and  the  calculation  may, 
thus,  be  abandoned  for  lattice  extensions  beyond  this  point. 

Subsequent  analysis,  in  pursuit  of  the  Griineisen  function,  will  also  require  the  derivative 
of  this  aggregated  stiffness  with  respect  to  lattice  spacing,  X.  Tedious,  but  straightforward, 
differentiation  gives 
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—(dFJdx)  .  (7-3) 
dk  T 


—OF  Jdx)  =  1/3  • 
dk  v,b 


dF  Jdx 

vib 


\5/2 


dFJdx 
\  L  J 


—  (dF./dx)  +  2/3  • 
dk  L 


d F  Jdx 

vib 


\5/2 


dFT/dx 


The  individual  L  and  T  component  and  aggregated  stiffness  values  may  be  obtained  via 
eqns  (5-6),  (6-5),  and  (7-2),  respectively.  In  addition,  derivatives  of  the  component  stiffnesses 
with  respect  to  k  are  obtained  by  taking  the  derivatives  of  eqns  (5-4)  and  (6-3)  with  respect  to 
k  (not  r). 

Though  the  angle,  y,  between  the  pairwise  force’s  line  of  action  and  the  direction  of 
vibration  does  not  vary  under  volumetric  change,  the  distance  sijk  does  change  according  to  the 
separation  metric  dsvoJdk,  eqn(4-3).  Thus,  to  accomplish  the  derivatives  of  the  stiffness 
components  required  of  eqn(7-3),  we  make  use  of  the  following  derivatives  of  eqns  (5-1)  and 
(5-2): 


d/dk(dsvib/dx)  =  0  ,  (7-4) 


and,  referring  to  eqn  (4-3), 


d/dk(d2svib/dxl)  =  -(sm2y)/5^  -dsvoltdk  .  (7-5) 

With  the  use  of  these,  the  derivatives  of  eqns  (5-4)  and  (6-3)  become 


d 

(  \ 
dFL 

00  00  00 

^  _  _  —  ujvo/ 

d2flds 2  \  cos2y  + 

r djlds  |s  < 
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dx 
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(i  *  0)  (7-6) 


and 
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sin2  y 

dk 

[dX  ) 

Z-j  r\ 

J  -  -oo  y-  -  oo  k~~  OO  Cl  A* 

S  ,  c2 

\  ) 

0  *0)  .  (7-7) 

It  may  again  be  confirmed  from  these  two  results  that,  were  the  lattice  composed  of  linear 
(harmonic)  springs  wherein  dflds  —  fls  and  d  2flds2  =  0,  the  variation  of  lattice  stiffness  with  X 
would  be  identically  zero.  These  two  equations  may  be  inserted  into  eqn  (7-3)  to  obtain  the 
derivative,  with  respect  to  X,  of  the  aggregated  stiffness.  Subsequent  use  may  be  made  of 
eqns  (7-2)  and  (7-3)  to  help  ascertain  the  vibrational  behavior  of  the  lattice.  The  availability  of 
aggregated  quantities,  however,  does  not  preclude  analysis  in  terms  of  the  individual  longitudinal 
and  transverse  ( L  and  T)  components  of  stiffness,  which  may  be  used  directly,  in  lieu  of 
eqns  (7-2)  and  (7-3),  respectively. 

8.  Results  and  Observations 


A  study  is  made  of  the  volumetric  and  vibrational  lattice  stiffnesses  using  the  preceeding 
equations.  The  lattice  structures  chosen  for  study  is  the  simple-cubic  (SC),  body-centered  cubic 
(BCC),  and  face-centered  cubic  (FCC)  structures.  The  sole  influence  of  the  BCC  and  FCC 
structures  on  the  equations  derived  above  is  that  z,  j,  and  k  necessarily  take  on  the  appropriate 
integer-plus-one-half  values  in  addition  to  integer  values.  The  method  described  here  could 
reasonably  be  extended  to  other  lattice  structures  as  well. 

To  proceed  with  the  analysis,  a  pairwise  potential,  e(s),  must  be  specified  with  energy 
expressible  as  a  function  of  separation  distance,  s,  so  as  to  determine  the  influence  of  any 
arbitrary  atom  (located  at  lattice  location  ijk)  upon  our  atom  of  interest,  O,  conveniently  located 
at  the  origin.  To  this  end,  we  select  the  so-called  H02  model  of  Schulte  and  Holzapfel  [15], 
primarily  because  they  noted  that  their  model,  in  contrast  to  others,  exhibited  the  proper  limiting 
behavior  in  strong  compression.  Though  this  model  was  not  originally  proposed  as  a  pairwise 
potential,  we  choose  it  nonetheless  for  study.  And  though  we  have  no  assurance  that  it  will  do - 
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so,  the  use  of  this  model  as  the  pairwise  potential  will  allow  it  to  be  determined,  not  only  if  the 
functional  form  of  the  aggregated  crystal  lattice  potential  follows  that  of  the  pairwise  interaction, 
but  whether  the  decay  rates  and  dimensions  of  the  many-atomed  crystal  lattice  are  similar  to 
those  of  the  corresponding  pairwise  interaction. 

Though  e(s)  is  taken  in  this  study  to  be  the  H02  model,  any  suitable  potential  may  be 
used  over  its  domain  of  compressive  accuracy.  The  far-field  behavior  of  the  potential  under  large 
expansions  is  significant  as  well  since  the  summations  of  the  lattice  quantities  are,  in  theory, 
taken  out  to  infinite  separation.  This  concern  is  especially  relevant  for  low  decay  rates, 
whereupon  the  potential  decays  slowly  with  increased  separation.  The  concern  is  not  one  of 
absolute  convergence,  since  the  number  of  atoms  in  the  summation  at  a  given  separation  distance 
goes  only  as  the  square  of  the  separation  distance,  whereas  the  potential  decays  exponentially. 
Rather,  the  concern  is  one  of  what  the  actual  converged  value  is;  a  cursory  examination  has 
shown  some  sensitivity  of  results  to  the  manner  in  which  the  potential  decays  with  separation 
( i.e .,  sensitivity  to  the  choice  of  pairwise  potential).  However,  this  sensitivity  has  been  one  of 
value,  not  of  trend. 

The  H02  potential,  being  used  to  describe  the  pairwise  interaction  of  any  two  atoms  in 
the  lattice,  actually  characterizes  the  force  interaction  instead  of  the  energy  potential.  Because 
the  H02  potential  is  being  used  to  characterize  the  pairwise  interaction,  rather  than  the  full  lattice 
interaction,  we  must  express  the  potential  in  terms  of  force  and  distance,  rather  than  pressure  and 
volume  (as  it  was  originally  given)  since,  for  a  pairwise  interaction,  pressure  and  volume  have 
no  meaning.  Since  F =pc  X2,  the  H02  expression,  expressed  in  terms  of  pairwise  force,  is 

/  =  s0(dflds)0  x'3  (1  -x)  exp  [a  (1  -x)]  ,  (8-1) 

where  a  is  the  interaction  decay  constant  of  the  pairwise  system  and  x  is  a  nondimensional 
interatomic  spacing  parameter,  given  by  x  =  s/s0.  The  parameter  s0(dflds)0  quantifies  the  force 
gradient  at  the  equilibrium  (zero-force),  pairwise  spacing,  s0.  The  pairwise  energy,  £(s),  and  the 
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force  derivative,  dfds,  may  be  directly  obtained  through  integration  and  differentiation  of 
eqn  (8-1),  respectively. 

This  study  is  performed  for  a  variety  of  representative  values  of  a,  with  results  being 
obtained  as  a  function  of  lattice  spacing  (i.e.,  the  governing  equations  are  summed  for  various, 
closely  spaced  increments  of  X/s0).  Since  summations  to  infinity  cannot  actually  be  performed, 
a  suitably  large  radius  of  influence,  smaj:,  is  chosen  for  the  summation  in  order  to  obtain 
convergence  ( i.e .,  terms  were  summed  if  For  typical  values  of  X  and  a,  values  of  smax 

equal  to  8  to  1 1  times  s0  were  observed  to  be  sufficient,  with  the  larger  values  corresponding 
to  smaller  values  of  a.  Of  course,  for  X  values  that  might  be  a  small  fraction  of  s0,  the  number 
of  summed  terms  grows  rapidly,  increasing  as  [(5^  /j0)(j0  /^)]3-  Furthermore,  as  the  pairwise 
decay  rate  is  successively  diminished,  such  that  the  effects  of  far-field  influence  decay  ever  more 
slowly  with  distance,  an  increasing  number  of  summation  terms  are  required  to  achieve 
convergence  at  a  single  highly  compressed  lattice  spacing. 


’  The  integration  of  the  H02  force  function,  eqn  (8-1),  by  standard  methods  requires  an  infinite  summation,  with  successive  terms  going  as 
x7( n  n\ ).  For  large  x ,  corresponding  to  large  lattice  spacings,  the  number  of  terms  required  for  convergence  and  the  limits  of  computer  precision 
make  this  approach  impractical.  Therefore,  a  Fit  was  developed  for  the  offending  exponential  integral.  By  making  use  of  the  substitution,  w  - 
l/(ax),  the  problem  may  be  expressed  in  terms  of 

-  I  <CU> 

’°  N"*'  • 

J  x  0 


The  fit  to  this  transformed  integral  has  been  found  to  be  accurate  to  within  0.008%  over  the  complete  domain  w  >  0  (0<  x  <  «=).  The  Fit  is  given 
by 


f 

w  -  In  1 1  * 

v 


-  [  w  -  ln(  1  ♦  w)  ] 


(1  ♦  3.918w  -  2.576w2) | 
(1  -  5.573 h-  5.888 w2y  j 


The  constants  in  this  Fit  are  tailored  to  minimize  the  relative  error  in  the  so-called  second  exponential  integral,  £;(x).  A  general  treatment  of  this 
approach  for  Fitting  the  exponential  integrals,  E„ ,  without  the  use  of  Fit  splicing,  may  be  found  in  Segletes  [27]. 


18 


A  number  of  general  observations  that  seemed  to  apply  regardless  of  the  value  of  a  were 
made  from  the  results  of  the  study  and  are  given  as  follows: 

(1)  The  equilibrium  condition  for  the  lattice  exhibits  a  nearest-neighbor  separation  distance 
less  than  that  for  the  pairwise  interaction.  In  particular,  the  equilibrium  spacing  of  the  nearest- 
neighbor  atoms  within  a  lattice,  seq,  must  always  be  smaller  than  the  pairwise  equilibrium 
spacing,  s0,  since  for  a  lattice  at  equilibrium,  the  attraction  of  the  far-field  atoms  must  be  exactly 
balanced  by  the  repulsion  of  the  nearest  neighbors.  These  nearest  neighbors,  on  the  other  hand, 
can  only  be  exhibiting  repulsion  if  their  separation  distance  is  smaller  than  the  equilibrium, 
pairwise  spacing.  Table  1  gives  the  equilibrium  lattice  separation  as  a  function  of  pairwise  decay 
rate  for  the  various  lattice  structures  (SC,  BCC,  and  FCC)  when  using  the  H02  potential  as  the 
pairwise  potential.  The  geometry  of  lattice  structures  other  than  simple  cubic  dictates  that  the 
interatomic  spacing  to  the  nearest  neighbor,  seq,  is  less  than  the  unit  cell  dimension,  A^,.  Thus, 
while  the  unit  cell  dimension  may  exceed  the  pairwise-equilibrium  spacing,  s0,  in  the  case  of 
BCC  and  FCC  structures,  the  nearest-neighbor  separation  distance,  seq,  will  not  do  so,  and  for  the 
aforementioned  reasons. 


Table  1.  Relative  Values  of  Equilibrium  Lattice  Spacing  as  a  Function  of  Pairwise 
Interaction  Decay  Rate 


a 

V*0  (Seq  Is o) 

SC 

BCC 

FCC 

i.i 

0.13  (0.13) 

0.16  (0.14) 

0.20  (0.14) 

1.5 

0.45  (0.45) 

0.55  (0.48) 

0.69  (0.49) 

2 

0.65  (0.65) 

0.79  (0.69) 

1.00  (0.71) 

2.5 

0.76  (0.76) 

0.92  (0.79) 

1.15  (0.82) 

3 

0.83  (0.83) 

0.99  (0.86) 

1.24  (0.88) 

4 

0.90  (0.90) 

1.06  (0.92) 

1.33  (0.94) 

5 

0.94  (0.94) 

1.09  (0.95) 

1.37  (0.97) 

oo 

1.00  (1.00) 

2/3 1/2  =  1.15  (1.00) 

2I/2  =  1.41  (1.00) 
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Also,  the  lower  the  decay  rate,  the  higher  the  equilibrium  compaction  (i.e.,  small  a 
implies  small  A^  /s0).  The  reason  for  this  behavior  is  that  a  smaller  exponential  decay  rate  in  the 
potential  implies  a  larger  number  of  far-field  neighbors  exerting  a  significant  attractive  force. 
This  increased  attractive  force  from  the  far  field  must  be  balanced  by  a  larger  repulsion  from  the 
nearest  neighbors.  This  far-field  attraction,  therefore,  causes  the  equilibrium  lattice  spacing, 
to  compress  significantly  with  respect  to  the  pairwise  equilibrium  distance,  s0.  For  these  low- 
decay-rate  cases,  the  form  of  the  pairwise  potential  is  exercised  at  both  the  repulsive  and 
compressive  limits  of  accuracy.  The  low  values  of  A^  /s0  drive  the  nearest  neighbor  into  an 
extremely  compressed  state,  while  at  the  same  time,  the  low  decay  rate  allows  for  significant 
pairwise  interactions  from  atoms  located  a  far  distance  away.  The  number  of  atoms  playing  a 
significant  role  in  the  equilibrium  calculation  will  continue  to  increase  as  the  interaction  decay 
constant  is  decreased,  with  the  computational  burden  increasing  correspondingly.  But  such 
behavior  is  only  a  limiting  case,  and,  in  general,  a  complete  and  converged  analysis  for  a  given 
value  of  a,  over  an  extensive  range  of  X,  requires  mere  seconds  of  computation  on  today’s 
generation  of  scientific  workstations  (e.g.,  SGI  Octane). 

(2)  The  summed  lattice  force,  F,  like  the  given  pairwise  force,/,  may  be  nominally  described 
over  a  moderate  range  of  lattice  spacings  by  the  H02  lattice  potential  of  Schulte  and 
Holzapfel  [15],  with  a  for  the  lattice  taken  identically  to  the  pairwise  interaction.  There  was 
no  reason  to  have  assumed,  a  priori ,  that  the  summed  potential  would  take  on  a  similar 
functional  form  to  the  pairwise  potential  composing  it.  Table  1,  for  example,  has  already 
indicated  that  the  absolute  minimum  of  the  summed  potential  occurs  at  a  different  lattice  spacing 
than  the  pairwise  potential.  Nonetheless,  it  was  observed  that,  up  to  moderate  compressions  (e.g., 
for  energies  near  and  below  the  binding  energy),  the  lattice  force  follows  a  similar  functional 
form  to  the  pairwise  force,  including  behavior  according  to  an  identical  decay  constant  (i.e.,  the 
value  for  a  characterizing  the  lattice  appears  to  be  identical  to  that  characterizing  the  pairwise 
interaction).  At  very  high  compressions  (and  particularly  for  larger  values  of  a),  the  scaled 
pairwise-  and  lattice-force  may  begin  to  diverge  slightly.  Figure  la  shows,  in  nondimensional 
terms,  the  pairwise-  and  lattice-force  for  an  FCC  lattice  with  a  =  2.5. 
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Figure  1.  The  pairwise  and  lattice  behavior  of  an  FCC  lattice  with  a  pairwise  interaction 
decay  rate  of  a  =  2.5:  (a)  nondimensional  force;  (b)  energy  potential,  depicting 
differing  values  for  pairwise-  and  lattice-stiffness,  T). 
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(3)  The  effective  nondimensional  stiffness  of  the  lattice,  t|,  is  different  from,  and  less  than, 
the  stiffness  associated  with  the  pairwise  interaction,  ri^.  Though  the  lattice  behavior  is 
nominally  describable  by  the  same  H02  function  as  the  pairwise  behavior,  the  one  notable 
difference  observed  at  equilibrium  between  the  pairwise  and  lattice  interaction  is  that  the  value 
of  the  lattice  stiffness,  h^{dFidk)0  differs  from  the  pairwise  stiffness,  s0(dfds)0.  One  good 
nondimensional  way  to  express  this  difference  is  by  way  of  the  i)  parameter  employed  by  Rose  et 
al.  [11].  For  the  pairwise  interaction,  it  is  given  by 

Tlpa.r  =  -idflds) o  s20/eh  ,  (8-2) 

where  eb  is  the  binding  energy  of  the  pairwise  interaction.  For  the  full  lattice,  the 
nondimensional  stiffness  is  given  by 

r|2  =  -3 (dF/dk)^  /  Eb  .  (8-3) 

In  eqn  (8-3),  Eb  is  the  binding  energy  of  the  lattice.  The  factor  of  3  reflects  the  definition  of 
lattice  force,  as  introduced  in  eqn  (4-4).  Both  the  pairwise  and  lattice  binding  energies  are 
obtainable,  from  the  integration  of  eqn  (8-1)  in  the  case  of  the  pairwise  interaction,  and  from  the 
triple  summation  of  eqn  (4—2)  for  the  lattice.  In  both  cases,  the  minimum  in  the  energy  potential 
is  the  binding  energy.  Table  2  shows  the  comparison  of  nondimensional  stiffnesses  (pairwise  and 
lattice)  arising  from  a  given  interaction  decay  rate  in  the  H02  potential.  As  a  is  lowered,  the  fact 
that  the  nondimensional  stiffness  of  the  lattice  decreases  relative  to  the  pairwise  stiffness  is 
indicative  of  the  fact  that,  for  low  values  of  a,  the  far-field  atoms  account  for  an  ever-increasing 
percentage  of  the  lattice’s  binding  energy. 

As  one  of  their  primary  contributions.  Rose  et  al.  [11]  did  note  the  functional  similarity 
between  the  energy  potentials  of  a  great  variety  of  bonding  configurations,  including  metallic  as 
well  as  diatomic.  And  though  the  pairwise  bond  in  a  crystal  lattice  configuration  is  not  identical 
to  the  corresponding  homopolar  bond  (because  of  differing  electronic  structures),  it  is  not  wholly 
unexpected  to  find  that  the  (summed)  lattice  potential  does,  indeed,  follow  a  form  similar  to  the 
pairwise  potential.  Figure  lb  depicts  the  lattice  potential  for  an  FCC  structure,  derived  via  the 
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triple  summation  of  eqn  (4—2),  using  a  pairwise  decay  rate  of  a  =  2.5,  corresponding  to 
r\pair  =  4.82  (see  Table  2).  Closely  tracking  that  curve,  shown  on  the  same  figure,  is  the 
(unsummed)  potential  of  Rose  et  al.  [11],  employing  a  value  of  rj  equal  to  4.55. 


Table  2.  Nondimensional  Stiffnesses3  as  a  Function  of  Interaction  Decay  Constant 


a 

pair 

n  (W>) 

sc 

BCC 

FCC 

i.i 

3.21 

2.65  (0.83) 

2.65  (0.82) 

2.65  (0.82) 

1.5 

3.70 

3.26  (0.88) 

3.25  (0.88) 

3.25  (0.88) 

2 

4.26 

3.94  (0.92) 

3.92  (0.92) 

3.92  (0.92) 

2.5 

4.82 

4.56  (0.95) 

4.55  (0.94) 

4.55  (0.94) 

3 

5.42 

5.15  (0.95) 

5.14  (0.95) 

5.14  (0.95) 

6.46 

6.27  (0.97) 

6.28  (0.97) 

6.28  (0.97) 

5 

7.45 

7.34  (0.98) 

7.37  (0.99) 

7.37  (0.99) 

3  The  pairwise  and  lattice  nondimensional  stiffnesses  were  acquired  using  eqns  (8-2)  and  (8-3),  respectively. 


One  point  of  note  regarding  eqn  (8-3)  is  that  the  parameter  t)  is  defined  in  terms  of  the 
slope  of  the  lattice  force  at  the  Xq  lattice-equilibrium  reference  state.  This  parameter  is  a 
nondimensional  measure,  therefore,  of  the  lattice  stiffness.  Nonetheless,  Rose  et  al.  [1 1]  referred 
to  the  parameter  t|  as  a  measure  of  the  lattice  anharmonicity.  Though  they  are  technically 
correct,  as  shall  be  shown,  we  prefer  here  instead  to  refer  to  T|  as  a  nondimensional  stiffness 
because  of  its  definition  in  eqn  (8-3). 

By  contrast,  the  term  anharmonicity  characterizes,  by  definition,  the  curvature  (and  not 
the  slope)  of  the  lattice  force  function  because  of  the  fact  that  the  force-curvature  for  a  harmonic 
spring  is  identically  zero.  If  we  understand  the  Griineisen  function,  T,  defined  in  eqn  (1-1),  to 
be  a  measure  of  the  vibrational  anharmonicity  of  the  lattice,  then  we  can  characterize  the 
volumetric  anharmonicity,  call  it  TV(,;,  in  terms  of  a  definition  similar  to  eqn  (1-1).  Relating  to 
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back  to  the  vibrational  stiffness  by  way  of  eqn  (2-2)  and  expressing  the  volume-derivative  in 
terms  of  A  instead,  eqn  (1-1)  may  be  expressed  as 

A3/y  =  r  =  -(A/6)  [d/dX(dFvlh  /a*)]  /  0 Fvlh  /dx)  .  (8-4) 

A  corresponding  definition  for  volumetric  anharmonicity,  employing  the  volumetric  rather  than 
the  vibrational  stiffness,  would  follow  as 

r,0/  =  -(A/6)  ( d2F/dX 2)  /  (dFIdX)  .  (8-5) 

Were  one  to  convert  A  back  into  V,  and  substitute  F  =  pcX2=pcV1},  eqn  (8-5)  could  be  shown  to 
be  precisely  the  relationship  derived  by  Dugdale  and  MacDonald  [5]  for  T,  who  had  assumed 
(along  with  many  others)  that  the  volumetric  and  vibrational  stiffnesses  were  functionally 
identical.  Table  3  provides  a  comparison  of  these  two  metrics  under  various  conditions.  Note, 
however,  that  T  in  the  present  study  is  composed  only  of  the  longitudinal  and  transverse 
stiffnesses  along  the  principal  directions  of  the  lattice,  whereas  a  full  accounting  of  T  would 
include  the  effects  from  all  lattice  modes.  We  know  from  the  work  of  Dugdale  and 
MacDonald  [5],  for  example,  that  measured  values  for  vibrational-  and  volumetric-r0,  for  real 
materials,  are  generally  closer  to  each  other  than  is  reflected  in  Table  3.  However,  though  the 
value  for  T0  may  be  incorrect  as  a  result  of  accounting  for  the  vibrational  modes  only  along  the 
principal  direction  of  the  lattice,  it  may  still  be  reasonable  to  assume  that  the  manner  in  which 
these  modes  functionally  change  with  lattice  spacing  is  characteristic  of  how  the  other  modes 
change. 

If  an  energy  potential,  one  that  precisely  satisfied  Rose’s  energy  potential,  were  evaluated 
in  eqn  (8-5)  at  the  reference  X0  state,  the  derived  value  of  rvoj0  would  be  q/2.6  (or  2.3r)/6,  to  be 
precise).  We  therefore  see  why  Rose  et  al.  referred  to  T)  as  an  anharmonicity  parameter,  even 
though  its  actual  definition,  eqn  (8-3),  characterizes  a  nondimensional  stiffness:  because  of  the 
particular  form  of  their  energy  potential,  the  parameter  r|  (or  more  accurately  t|/2.6)  characterizes 
the  volumetric  anharmonicity,  as  well. 
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Table  3.  Equilibrium  Values  for  the  Vibrational3  and  Volumetric  T 


a 

r0  (rvo;0 ) 

scb 

BCC 

L  . . 

FCC 

i.i 

0.72  (1.37) 

0.70  (1.37) 

0.70  (1.36) 

1.5 

0.89  (1.49) 

0.83  (1.48) 

0.83  (1.48) 

2 

1.10  (1.65) 

1.01  (1.64) 

1.00  (1.64) 

2.5 

1.31  (1.82) 

1.19  (1.81) 

1.17  (1.81) 

3 

1.52  (1.98) 

1.38  (1.97) 

1.36  (1.97) 

4 

1.92  (2.32) 

1.74  (2.31) 

1.72  (2.31) 

5 

2.32  (2.65) 

2.09  (2.65) 

2.09  (2.65) 

Vibrational  T  computed  from  principal  stiffnesses  only. 

b  Vibrational  F  for  simple  cubic  lattice  computed  from  longitudinal  stiffness  only,  as  SC  lattice  not  stable  in  shear,  along 
the  principal  lattice  directions. 


Even  though  the  triple-summed  energy  potential  in  the  current  effort,  eqn(4-2),  is 
nominally  describable  by  Rose’s  potential  (Figure  lb),  it  would  be  presumptious  to  conclude  that 
rvo/0  will  always  equal  rj/2.6.  Therefore,  we  should  use  TvM  rather  than  tj  to  characterize  the 
volumetric  anharmonicity.  For  the  case  of  pairwise  interactions  only,  where  the  H02  potential 
has  been  assumed  as  the  governing  form,  the  relation  rvo;0  =  1  +  a/3  governs  the  pairwise 
curvature  at  the  equilibrium  state. 

The  volumetric  anharmonicity,  rvo/0,  since  it  defines  the  curvature  of  the  lattice  force  at 
the  Xq  state,  should  also  be  relatable  to  the  parameter  B',  which  is  employed  by  a  number  of 
authors  [15-19],  and  represents  the  pressure  derivative  of  the  bulk  modulus,  B,  as 

B1  =  (BB/dp)T  ,  (8-6) 

though  some  authors  have  defined  this  derivative  as  being  at  constant  pressure,  not  temperature. 
The  confusion  arises  because  eqn(8-6)  defines  a  curvature  along  the  reference  isotherm. 
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However,  actual  data  tend  to  be  obtained  from  a  number  of  data  points  at  constant  pressure.  The 
expression  B'0  refers  to  this  quantity  at  the  reference  state  of  zero  temperature  and  pressure.  If 
we  expand  the  definition  of  B  as  -V  ( dp/dV)T ,  then  eqn  (8-6)  may  be  expressed  as 

b'  =  -i  -  voyavV(dp/avor  •  (8-7) 

Eqn  (8-7)  looks  very  similar  to  the  definition  of  T  proposed  by  Slater  [2],  TSh  and  may  be  related 
to  it  as 

r5,  =  -1/6  +  572  .  (8-8) 

A  relation  similar  to  this  was  actually  posited  by  Syassen  and  Holzapfel  [  1 8],  though  their 
relation  was  based  on  model  calculations,  and  not  derived.  At  the  reference  state  only,  the  value 
for  r0  predicted  by  Slater’s  model  is  1/3  in  excess  of  the  value  later  predicted  by  Dugdale  and 
MacDonald’s  model  [5],  TDW0  =rs(0- 1/3.  And  since  eqn  (8-5)  indicated  that  the  Dugdale- 
MacDonald  TDM  was,  in  fact,  Tvol,  we  therefore  have  that 

rvol0  =  -1/2  +  B'Jl  .  (8-9) 

This  equation  constitutes  yet  one  more  method  for  estimating  Tl(>(0,  though  there  appears  to  be 
wide  scatter  in  the  experimental  measurement  of  B'0,  more  so  than  in  the  estimation  of  T|. 

(4)  The  vibrational  and  volumetric  spring  constants  are  quantitatively  distinct  functions. 

Perhaps  a  primary  result  of  this  analysis  is  quantitatively  verifying  the  distinct  nature  of  the 
volumetric  and  vibrational  stiffness  of  an  atomic  lattice,  a  distinction  which  has  been  recently 
asserted  on  qualitative  grounds  [13]  but  one  which  contradicts  the  assumption  employed  by  many 
previous  analyses  on  the  subject  [2, 5, 6-9].  Since  the  distinction  is  shown  to  arise  from  the 
influence  of  nonnearest  neighbors,  we  might  expect  the  magnitude  of  the  discrepancy  to  be 
greater  for  low-decay-rate  interactions.  Figures  2  and  3  bear  this  hypothesis  out  by  depicting  a 
comparison  of  the  volumetric  and  vibrational  stiffnesses  for  two  FCC  lattices  of  differing 
stiffnesses.  In  terms  of  absolute,  not  relative,  quantities,  Figure  4  displays  the  energy,  force,  and 
stiffnesses  of  an  FCC,  a  =  5  lattice. 
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Figure  2.  A  relative  comparison  of  the  volumetric  stiffness,  dF/dk,  to  the  vibrational 
stiffness,  dFvii/dx,  for  an  FCC  lattice  with  a  pairwise,  interaction  decay  rate  of 
a  =1.5:  (a)  to  2.5  times  the  nominal  stiffness;  (b)  to  10  times  the  nominal 

stiffness. 
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Figure  3.  A  relative  comparison  of  the  volumetric  stiffness,  dF/dk,  to  the  vibrational 
stiffness,  dFvli  fdx,  for  an  FCC  lattice  with  a  pairwise,  interaction  decay  rate  of 
a  =  4. 


XA0 


Figure  4.  The  energy,  force,  and  volumetric  and  vibrational  stiffnesses  for  an  FCC  a  =  2.5 
lattice,  shown  in  comparison  to  each  other  (note  that  the  stiffness  functions  are 
shown  as  positive  in  compression  for  convenience  in  visualization). 
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(5)  The  zero  of  the  vibrational  spring  constant  is  observed  to  be  at  a  value  of  X  greater  than 
the  zero  of  the  volumetric  spring  constant.  As  a  decreases,  the  distance  between  these 
zeroes  increases  as  well.  For  low  values  of  pairwise  stiffness,  the  zero  of  the  vibrational 
stiffness  can  be  significantly  different  from  that  of  the  volumetric  stiffness.  For  the  a  =  1.5  case 
depicted  in  Figure  2,  the  zeroes  of  the  volumetric  and  vibrational  stiffness  occur  at  X/Xq  values 
of  1 .28  and  2.41 ,  respectively.  This  difference  in  lattice  spacing  corresponds  to  a  huge  difference 
in  relative  volumes,  V/V0,  of  2.1  and  14.0,  respectively.  At  a  large  enough  expansion,  however, 
even  the  vibrational  stiffness  is  observed  to  become  negative.  When  the  vibrational  stiffness 
becomes  negative,  the  characteristic  frequency,  (0,  becomes  undefined  and  the  Griineisen 
relationship,  eqn(l-l),  breaks  down.  The  zero  of  the  volumetric  spring  constant  (i.e.,  the 
inflection  point  of  the  energy  potential)  has  been  associated  with  a  metal’s  melting  point  [9],  and 
perhaps  it  is  this  phase  change  to  a  liquid,  not  considered  in  this  study,  that  in  reality  precludes 
the  Griineisen  singularity  associated  with  the  zero  of  the  vibrational  spring  constant. 

(6)  The  value  of  the  Griineisen  function,  T,  is  observed  to  approach  2/3  in  the  high- 
compression  limit.  By  contrast,  r„o;  does  not.  Eqn  (8-4)  shows  how  the  Griineisen  function 
(and  x\r  =  WT)  may  be  evaluated  when  the  vibrational  stiffness  is  known.  The  stiffness  itself  is 
obtained  from  eqn  (7-2)  and  its  derivative  with  respect  to  X  from  eqn  (7-3).  Alternately, 
longitudinal-  and  transverse-component  values  of  T  may  be  obtained  by  using  the  L  and/or  T 
component  values  directly  in  eqn  (8-4).  The  value  for  \jr  at  any  given  value  of  X  is  WT, 
alternately  X3/T.  A  typical  example  of  the  T  function,  as  computed  in  this  maimer,  is  shown  in 
Figure  5  (expressed  in  X,  not  V),  along  with  \|r.  For  illustration,  both  the  L  and  T  component 
contributions  are  shown  as  well. 

The  functional  form  for  \|/  has  been  observed  to  be  similar  in  appearance  (to  that  of 
Figure  5)  over  the  full  range  of  cases  tested.  Therefore,  we  seek  to  further  examine  the  nature 
of  this  function.  The  H02  potential,  as  noted  by  the  authors  [15],  exhibits  the  proper  functional 
form  in  the  high-compression  limit.  We  would,  therefore,  hope  that  the  evaluation  of  \|f  and  T 
by  way  of  eqn  (8-4)  would  do  the  same.  In  fact,  we  note  here  that,  regardless  of  the  value 
employed  for  a,  or  the  lattice  structure  employed,  the  limiting  high-compression  value  for  T 
approaches  the  recognizable  value  of  2/3.  By  contrast,  were  there  no  distinction  made  between 
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Figure  5.  The  T  and  corresponding  \jr  functions  for  an  a  =  2.5,  BCC  lattice,  showing  both 
the  L  and  T  components,  as  well  as  the  aggregated  function. 


30 


the  vibrational  and  volumetric  stiffness,  then  rv£>/,  evaluated  by  way  of  eqn  (8-5),  should  instead 
provide  the  proper  limiting  behavior  of  2/3.  However,  rv0/  has  been  observed  not  to  approach 
this  limiting  value.  The  only  way  that  Tto/  could  approach  the  proper  limiting  value,  given  the 
H02  pairwise  potential,  would  be  if  interactions  beyond  the  nearest  neighbors  were  ignored. 
Unfortunately,  no  matter  how  large  an  interaction  decay  constant,  a,  might  be  chosen  in  an 
attempt  to  eliminate  the  influence  of  all  but  the  nearest  of  neighbors,  the  nearest  nonneighbors 
will  nonetheless  have  an  effect  in  the  high-compression  limit,  since  no  matter  how  much  a  is 
increased,  X  may  shrink  toward  zero  a  corresponding  amount  in  order  to  achieve  the  same 
relative  influence  of  its  nonnearest  neighbors. 

Considering  the  other  end  of  the  domain,  at  large  X,  a  negative  slope  on  \|/  has  been 
associated  with  one  type  of  thermodynamic  instability  at  elevated  temperatures  [20].  We, 
therefore,  examine  the  state  of  the  lattice  interactions  at  this  critical  point  where  dyldV = 0.  If 
we  denote,  at  this  maximum  in  the  \|r  curve,  the  lattice  spacing  as  Xslab  and  the  distance  to  an 
atom’s  nearest  neighbor  in  the  lattice  as  sslab,  then  we  may  learn  something  by  the  behavior  of 
these  variables  as  a  function  of  pairwise  stiffness  and  lattice  structure.  Table  4  portrays  these 
variables  in  a  normalized  form.  The  results  show  behavior  that  is  somewhat  uniform  across 
lattice  structures  (allowing  for  the  fact  that  the  SC  lattice  calculations  are  based  on  longitudinal 
stiffness  only).  Examining  sstab,  the  values  relative  to  s0  are,  in  addition,  nearly  independent  of 
the  interaction  decay  rate,  a.  One  may  conclude  from  this,  regardless  of  lattice  structure  and 
stiffness,  that  the  extremum  in  the  \|r  curve  occurs  when  the  interatomic  distance  to  the  nearest 
neighbor  in  the  lattice  is  approximately  equal  to  the  equilibrium  spacing  associated  with  the 
pairwise  potential.  At  such  a  lattice  spacing,  this  nearest  neighbor  is  barely  applying  any  force 
at  all,  and  all  neighbors  beyond  s0  are  applying  an  attractive  force.  Thus,  the  maximum  of  the 
\jf  function  occurs,  for  a  variety  of  lattice  structures  and  over  a  large  range  of  nondimensional 
stiffness,  at  nearly  the  identical  relative  pairwise  configuration — a  configuration  associated  with 
a  minimal  pairwise  interaction  from  the  nearest  neighbor. 

Considering,  however,  the  situation  in  terms  of  the  lattice  behavior,  as  opposed  to  the 
pairwise  interactions,  an  examination  of  Xs!ab  values  from  Table  4  reveals  where  the  equilibrium 
spacing  for  the  lattice  falls  relative  to  the  critical  spacing.  At  low  decay-rate  values,  the  critical 
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value  for  X  falls  at  several  times  the  equilibrium  spacing  for  the  lattice.  As  the  decay  rate  is 
increased,  however,  Table  4  reveals  that  the  equilibrium  lattice  spacing  approaches  the  critical 
value  for  X. 


Table  4.  Critical  Lattice  Spacings  Associated  with  dy/dV  =  0a 


a 

Xitab  ( Sstab  /S0) 

scb 

BCC 

FCC 

i.i 

7.09  (0.90) 

7.34(1.00) 

7.39  (1.03) 

1.5 

2.04  (0.92) 

2.10  (1.00) 

2.12  (1.04) 

2 

1.42  (0.93) 

1.46  (1.00) 

1.48  (1.04) 

2.5 

1.22  (0.93) 

1.26  (1.00) 

1.27  (1.04) 

3 

1.13  (0.93) 

1.16  (1.00) 

1.17  (1.03) 

4 

1.04  (0.93) 

1.07  (0.98) 

1.07  (1.01) 

5 

1.00  (0.93) 

1.02  (0.97) 

1.02  (1.00) 

a  The  v  function  computed  from  principal  stiffnesses  only. 

b  The  \]/  function  for  simple  cubic  lattice  computed  from  longitudinal  stiffness  only,  as  SC  lattice  is  not  stable  in  shear, 
along  the  principal  lattice  directions. 
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PART  II — Macroscopic  Lattice  Behavior 


9.  Mechanical  Characteristics  of  the  Lattice 


The  summed  lattice  equations  that  have  been  derived  in  Part  I  [eqns  (4-2),  (4-7),  (4-8), 
(7-2)  and  (7-3)]  depict  a  methodology  for  modeling  the  interatomic  forces  and  vibrational 
characteristics  of  an  atomic  lattice.  As  a  paradigm,  it  is  simple,  complete,  and  of  general 
applicability.  Unfortunately,  from  a  practical  point  of  view,  a  triple  summation  is  an  inefficent 
vehicle  by  which  to  computationally  describe  the  behavior  of  a  lattice.  Furthermore,  the  insights 
provided  by  the  triple  summation  are  limited  in  scope. 

It  was,  therefore,  a  goal  in  this  study  to  capture  the  essence  of  the  governing  equations, 
which  involve  triple  summations,  in  terms  of  a  simpler  analytical  formulation.  Our  approach  is 
to  examine  separately  the  mechanical  and  the  thermal  characteristics  of  the  summed  lattice 
equations,  each  in  terms  of  the  lattice’s  vibrational  behavior.  These  separate  mechanical  and 
thermal  models  can  then  be  combined  into  a  unified  equation  of  state. 

Based  in  part  on  the  semi-empirical  form  of  eqn  (2-1),  similar  looking  forms  were  tested 
in  an  attempt  to  analytically  describe  the  summed  lattice  force  without  the  use  of  fitting 
parameters.  A  remarkably  good  fit  was  found  possible,  over  a  wide  range  of  stiffnesses  and 
lattice  spacings,  with  the  following  form: 


F 

-\(dF/dk)0 


- (cd/o)Jp  ln((o/con)  . 

3T0(^0)2 


(9-1) 


This  form  is  guaranteed  to  provide  slope  compatibility  of  the  force  at  Xq,  while  the  parameter  P 
is  selected  to  guarantee  curvature  compatibility  in  the  lattice  force  at  Xq.  That  requirement 
dictates  that 
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P  =  [r -  1/3  -  1/2  -r0(<fy/<Ao0]/r0  ♦ 


(9-2) 


Furthermore,  we  note  that  3  is  obtainable  from  material  properties  at  the  force-free  equilibrium 
state. 


Let  us  recall  that  given  in  this  problem  is  a  known  pairwise  interaction  potential, 
integrated  from  eqn  (8-1),  which  is  characterized  by  the  parameters  a,  s0,  and  s0(d/Ids)Q.  Through 
direct  differentiation  of  this  pairwise  force,  the  pairwise  stiffness  may  be  readily  obtained.  In 
eqn  (9-1),  which  we  propose  to  employ  to  characterize  the  behavior  of  the  full  lattice,  and 
XoidF/dX) o  are  computed  from  the  equilibrium  condition  of  eqn  (4-2).  The  parameter  3  may  be 
computed  from  eqn  (9-2).  The  ratio  (co/co0),  from  the  definition  in  eqn  (2-2),  is  given  as 

to/CDo  =  [dF^  /dx  /  0 Fvlh  /a*)0  ] 1 2  .  (9-3) 

Both  of  the  constituent  terms  composing  the  right-hand  side  of  eqn  (9-3)  are  obtained  by  way 
of  eqn  (7-2).  The  value  for  T0  may  be  obtained  from  eqn  (8-4),  being  evaluated  at  the 
equilibrium  value  of  lattice  spacing. 

Thus,  we  are  able  to  evaluate  eqn  (9-1)  without  the  use  of  fitting  parameters  and  compare 
the  result  to  the  force,  F,  obtained  by  way  of  the  triple  summation  of  eqn  (4-7).  Figure  6  shows 
examples  of  the  quality  with  which  eqn  (9-1)  may  describe  the  lattice  force  over  a  wide  range 
of  interatomic  decay  rates.  The  fits  are  for  an  FCC  crystal,  though  the  other  crystal  types  (e.g., 
SC  and  BCC)  show  identical  trends  with  comparable  accuracy.  The  comparisons  are  shown  to 
a  compression  corresponding  to  a  stiffness  of  15x  the  equilibrium  stiffness,  which,  in  all  cases 
shown,  corresponds  to  a  maximum  relative-volumetric  compression  at  or  beyond  V7  V0  =  0.22 
(AAfl  at  or  beyond  0.6).  The  comparison  is  of  excellent  quality  for  the  low-  and  mid-range  decay 
rates,  though  it  does  begin  to  deviate  somewhat  for  the  a  =  5  case. 

It  should  be  noted  that,  unlike  eqn  (4-7),  which  is  valid  for  all  A  states,  the  analytical, 
frequency-based  eqn  (9-1)  can  only  apply  where  a  vibrational  frequency,  to,  is  definable.  When 
dFvib  /dx  changes  sign  (at  some  expansion  beyond  the  inflection  point  in  the  lattice  potential, 
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Figure  6.  A  comparison  of  the  numerically  computed  (summed)  lattice  force  [eqn  (4-7), 
solid  line]  to  that  analytically  evaluated  by  way  of  eqn  (9-1)  (dashed  line),  for  an 
FCC  lattice  of:  (a)  a=  1.1;  (b)  a=  3;  (c)  a=  5. 
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(C) 

Figure  6.  A  comparison  of  the  numerically  computed  (summed)  lattice  force  [eqn  (4-7), 
solid  line]  to  that  analytically  evaluated  by  way  of  eqn  (9-1)  (dashed  line),  for  an 
FCC  lattice  of:  (a)  a  =  1.1;  (b)  a  =  3;  (c)  a  =  5. 
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eqn  (9-1)  can  no  longer  apply.  Of  course,  by  this  level  of  expansion,  dF/dk  has  already  changed 
sign,  whereupon  the  lattice  can  no  longer  remain  mechanically  stable  ( i.e .,  any  slight  deviation 
from  the  ideally  spaced  lattice  configuration  will  result  in  negative  restoring  forces,  which  will 
effectively  tear  apart  the  lattice). 

10.  Thermal  Characteristics  of  the  Lattice 

Whereas  the  interatomic  potential  is  responsible  for  the  lattice’s  zero-temperature,  force 
of  repulsion,  as  computed  by  eqn  (4-7),  thermal  excitations  also  add  to  the  net  interatomic 
repulsion.  Though  temperature  is  a  direct  measure  of  the  thermal  excitations,  the  internal  energy 
of  thermal  excitation,  Eth  -  (is  -  Epot)v,  is  directly  related  to  the  temperature  by  way  of  the  specific 
heat  function  and  is,  therefore,  another  measure  of  the  thermal  excitation.  As  noted  in  this 
paper’s  introduction,  it  is  the  Gruneisen  function,  T,  that  directly  relates  a  material’s  thermal 
energy  change  to  a  change  in  pressure.  Knowledge,  therefore,  of  the  T  (or,  alternately,  t|f) 
function,  in  addition  to  the  lattice  potential,  constitutes  a  complete,  thermal  equation  of  state. 
Note,  however,  that  because  the  vibrational  stiffness  modeled  in  this  paper  is  not  considered  to 
be  a  function  of  vibrational  amplitude  [i.e.,  infinitesimal  amplitudes  have  been  assumed  as  part 
of  eqn  (2-2)],  the  current  analysis  will  necessarily  compute  T  as  a  function  of  volume  only. 
Fortunately,  such  an  assumption  is  widely  accepted  in  the  community,  over  a  large  range  of 
temperatures. 

In  addition  to  being  analytically  related  to  co  by  way  of  eqn  (1-1),  it  was  noted  that  the 
function  could  be  fit  to  to  in  another  manner,  again  over  a  wide  range  of  lattice  spacings  and 
interaction  decay  rates.  The  form  of  the  fit  is  given  by 

\|r/A4  =  iir/V*3  ~  w1  ,  (10-1) 

where  ^  is  a  parameter  that  can  be  fitted  for  a  given  lattice  configuration  and  decay  rate. 
Figure  7  provides  an  example  of  the  quality  of  this  correlation.  Similarly  good  correlations  were 
noted  for  other  lattice  types  and  decay  rates,  with  values  for  approaching  0.50  for  low 
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stiffnesses  (a~  1.1)  and  rising  to  approximately  0.56  for  high  stiffnesses,  regardless  of  lattice 
structure  (a  ~  5). 

Though  we  know  of  no  theoretical  reason  why  the  correlation  described  by  eqn(lO-l) 
should  hold  true,  a  further  measure  of  the  quality  of  the  correlation  may  be  obtained  by 
differentiating  eqn(lO-l)  with  respect  to  V  and  dividing  the  result  by  eqn(lO-l)  itself. 
Substitution  of  eqn  (1-1)  into  the  result  allows  for  the  elimination  of  to,  to  obtain 

£  =  4/3  -y/V-dy/dV  .  (10-2) 

Effectively,  eqn  (10-2)  compares  not  the  values  of  the  terms  composing  eqn  (10-1),  but  rather 
their  derivatives,  and  is  thus  a  more  stringent  comparison  of  the  correlation.  Figure  8  depicts, 
for  the  same  lattice  configuration  employed  in  Figure  7,  the  two  terms  on  the  right-hand  side  of 
eqn  (10-2).  Based  on  the  fitting  form,  the  vertical  distance  between  the  two  curves  should  be 
a  constant,  and  is  a  measure  of  The  correlation  is  seen  to  hold  over  the  compressive  domain 
and  only  begins  to  diverge  at  the  larger  lattice  spacings,  near  and  beyond  the  critical  value  of  X, 
where  dy/dk  changes  sign. 

Furthermore,  by  evaluating  and  its  derivative  at  the  equilibrium  lattice  spacing,  Xq, 
eqn  (10-2)  provides  a  nonempirical  way  to  select  a  value  for  Of  course,  evaluating  ^  in  this 
manner  does  not  necessarily  provide  the  best-fit  for  \\r  over  a  wide  range  of  X — however,  it  does 
guarantee  that  the  analytical  function  represented  by  eqn  (10-2)  will  match  the  value  and  slope 
of  \|/  at  Xq,  otherwise  obtained  from  the  triple-summation  method.  Whereas  a  good  fit  to 
eqn  (10-1)  over  a  wide  range  of  X  was  obtained  in  Figure  7  with  a  value  of  ^  =  0.50,  the  use  of 
eqn  (10-2)  to  analytically  evaluate  |  at  the  X^  reference  condition  provides  a  slightly  altered  value 
of  |  =  0.538,  which  we  use  for  the  remainder  of  this  comparison. 

A  solution  to  the  differential  equation,  eqn  (10-2),  is 

y/V„  =  (V/  Vo)  -  (3|  -  l/r0)  (V/  V</3  .  (10-3) 
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When  we  compare  this  explicit  form  for  to  that  computed  from  the  pairwise  summation  by 
way  of  eqns  (7-2),  (7-3),  and  (8-4),  we  fmd,  in  Figure  9a,  a  good  match  into  the  compressive 
regime,  with  a  divergence  noted  in  the  expanded  states  of  the  lattice,  near  and  beyond  the 
maximum  in  x\r.  If  we  invert  \\r  from  Figure  9a  and  express  the  result  in  terms  of  T,  as  in 
Figure  9b,  we  see  that  at  high  compressions,  results  may  deviate  as  well.  The  reason  here  is 
clear:  in  order  that  =  2/3,  must  necessarily  take  on  a  value  of  1/2.  As  noted  previously, 
values  for  required  to  fit  eqn(lO-l)  fall  in  the  0.5  to  0.56  range.  What  this  says  is  that 
eqn  (10-1)  is,  at  best,  a  fit  but  that  the  range  of  lattice  spacing  over  which  this  fit  applies  can  be 
quite  large,  indeed.  Recall,  as  well,  that  V/V0=(X/X0)3,  and  so  divergence  of  the  fit  at  a  relatively 
small  value  of  X  corresponds  to  an  extremely  compressed  state  in  terms  of  V. 

Eqn  (1-1)  may  be  integrated,  in  light  of  eqn  (10-3),  to  recover  a  functional  form  for  the 
lattice  frequency.  Again,  like  eqns  (10-2)  and  (10-3),  this  computed  frequency  is  based  upon 
the  correlation  cited  in  eqn  (10-1).  The  result  is 

co/co0  =  [3T0£(V^-1)  +  l]171  .  (10-4) 

Figure  9c  shows  the  analytically  computed  frequency  of  eqn  (10-4)  compared  to  the  summed 
value  for  (o/(D0,  acquired  from  eqn  (9-3)  by  way  of  the  triple  summation  of  eqn  (7-2).  To  more 
fully  ascertain  the  effect  of  the  divergence  noted  in  Figures  8  and  9a,  we  zoom  in  on  the  low- 
frequency,  expanded  lattice  region  of  Figure  9c  and  present  the  result  as  Figure  9d. 

Even  though  the  parameter  ^  may  be  obtained  from  known  material  properties  at  the  Xq 
reference  state,  the  fact  remains  that  eqn  (10-3)  is  still  a  simplified  fit  to  the  actual  functional 
form  exhibited  by  the  \|/  variable.  At  large  compressions,  the  form  of  eqn  (10-3)  may  diverge 
from  data.  At  large  expansions,  as  well,  there  is  a  rather  notable  departure  of  eqn  (10-3)  from 
the  actual  value  for  \|f.  Because  values  for  Xstab  l\,  as  given  in  Table  4,  are  large  for  small  a, 
the  analytical  equations  deriving  from  the  fit  to  eqn  (10-1)  should  be  accurate  at  the  nominal 
lattice  spacing  and  to  quite  an  extent  into  the  expanded  regime  for  small  a.  By  contrast,  for 
large  values  of  a,  the  analytical  fit  associated  with  eqn  (10-3)  will  be  valid  over  a  much  smaller 
range  of  expanded  states  because  the  Xq  reference  state  is  at  the  point  of  greatest  curvature  in  the 
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Figure  9.  A  comparison  of  the  computed  (summed)  lattice  properties  to  the  analytical  tit 
implied  by  eqns(10-3)  and  (19-4)  for  an  a  =  2.5,  BCC  lattice:  (a)  v/V0;  (b)  T; 
(c)  (o/(O0  (compression);  and  (d)  u/(O0  (expansion). 


(c) 


(d) 

Figure  9.  A  comparison  of  the  computed  (summed)  lattice  properties  to  the  analytical  fit 
implied  by  eqns(10-3)  and  (10-4)  for  an  a  =  2.5,  BCC  lattice:  (a)  \| //V0;  (b)  T; 
(c)  cd/cd0  (compression);  and  (d)  co/o)0  (expansion). 
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\|/(X)  curve.  Clearly,  the  eventual  solution  to  this  dilemma  is  a  better-fitting  functional  form  for 
\|/,  in  preference  to  eqn(10-3).  Nonetheless,  the  simplicity  offered  by  eqn(lO-l)  and  the 
equations  derived  from  it  will  help  us  to  formulate,  in  the  next  section,  a  complete,  thermal 
equation  of  state,  capable  of  modeling  general  thermodynamic  transitions  and  not  only  zero- 
temperature  ones. 

Interestingly,  in  the  low-a  limit,  where  £,  approaches  1/2  and  T0  (see  Table  3)  approaches 
2/3,  we  have  from  eqns  (10-3)  and  (10-4)  that  3/2  •  V  and  w  -  X"2.  The  conclusions,  in  this 
limiting  case  only,  are  that:  T  is  a  constant,  equal  to  2/3;  Xs!ah  approaches  infinity;  and  the 
vibrational  stiffness  never  changes  sign,  such  that  to  is  defined  out  to  infinite  lattice  separation. 

11.  A  Frequency-Based  Analytical  Equation  of  State 


The  lattice  force  description,  given  by  eqn  (9-1),  has  been  shown  to  describe  the  behavior 
of  SC,  BCC,  and  FCC  lattices  over  a  wide  range  of  lattice  spacings  and  interatomic,  interaction- 
decay  rates,  a.  When  given  a  pairwise  potential,  no  fitting  parameters  are  required  in  order  to 
describe  the  mechanical  behavior  of  the  lattice.  By  contrast,  the  semi-empirical  equation-of-state 
form  given  by  eqn  (2-1),  similar  in  form  to  eqn  (9-1),  requires  one  fitting  parameter  to  describe 
the  functional  behavior  of  \\r,  since  there  is  no  pairwise  potential  from  which  eqn  (2-1)  derives. 
The  semi-empirical  form,  eqn  (2-1),  does,  however,  have  one  advantage  over  eqn  (9-1) — that  it 
may  be  readily  integrated,  because  of  the  \jr  term,  to  yield  the  lattice  potential  as 

Epo,  =  ~Eb  ((D/a),/  [1  -  ln(o)/co0)A  ]  ,  (11-1) 

whereas  eqn  (9-1)  cannot  be  integrated  without  first  specifying  the  manner  in  which  to  is  coupled 
to  X. 


However,  begining  with  the  lattice  force  given  by  eqn  (9-1)  and  substituting  the  observed 
correlation  given  by  eqn  (10-1)  into  it  [in  effect,  multiplying  by  unity  in  the  form  of 
(co/g)0)^  '(X/Xq)4/  (v|rAjr0)],  one  may  obtain 
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F 


(11-2) 


-\{dF!dk\ 


V(3ro)  •  *•’/¥  ’  (w/w0)^ln(w/w0)  . 


In  deriving  eqn(ll-2),  use  has  been  made  of  the  definition,  given  in  (1-1),  that  \|t  =  V/T.  A 
comparison  of  this  result  to  eqn(2-l)  reveals  identical  forms  with  but  one  exception — in  the 
semi-empirical  form,  the  exponent  on  the  lattice  frequency  term  is  K,  whereas  in  the  present 
analysis,  it  is  (^+  P).  Based  on  the  the  values  of  ^  and  P  determined  for  the  lattices  studied,  one 
would  expect  to  find  values  for  the  exponent  roughly  in  the  4/3  to  5/3  range.  By  contrast,  values 
for  K  display  a  range  of  2/3  to  4/3. 

The  reason  that  the  semi -empirical  eqn  (2-1)  retained  K  as  the  exponent  on  to  arose  from 
the  mistaken  belief  that  to  remained  definable  out  to  infinite  expansions  and  that  the  form  of  the 
potential,  eqn  (1 1-1),  should  continue  to  fit  the  data  as  X  °°.  In  contrast,  the  current  analysis, 
leading  to  the  explicit,  triple-summation  calculation  of  dFvib/dx  and  (0,  has  shown  that  at  large 
enough  expansions  the  vibrational  stiffness  changes  sign  and,  thus,  the  characteristic  lattice 
frequency  becomes  undefined  beyond  that  point.  Acknowledging,  therefore,  the  inability  of  any 
frequency-based  equation  of  state  for  solids  to  capture  the  highly  expanded  behavior  of  the  lattice 
potential,  the  requirement  for  retaining  the  exponent  on  to  equal  to  the  value  K  is  no  longer 
imperative. 

Unlike  eqn  (9-1),  eqn  (11-2)  now  becomes  directly  integrable  to  acquire  the  energy 
potential.  If  we:  (1)  denote  the  sum,  %  +  P,  as  k;  (2)  recall  that  F  =  pcX2;  (3)  express  -X^(dF/dX)Q 
in  terms  of  the  equilibrium  bulk  modulus,  B0,  as  -\{dFldX\  =  3 B0  ;  and  (4)  reference  the 
energy  zero  to  the  equilibrium  lattice  state,  such  that  Ec  =  Epot  +Eb;  then  the  lattice  pressure  and 
energy  may  be  expressed  as 


pc  =  FIX2  =  B0  V0  /Tl  ■  (cD/co0r  ln(co/(O0) /\|f  (11-3) 

Ec  =  B0  V0  /(r0  k)2  <1  -  (w/co0)K  [1  -  ln(&)/to0)K  ]  >  .  (11-4) 
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The  quantity  B0  V0,  appearing  only  as  a  product,  may  be  replaced  by  the  square  of  the 
equilibrium,  bulk  sound  speed,  C\.  Note  also  that,  unlike  eqn  (9-1),  the  H02  potential,  or  many 
of  the  other  popular  potentials  in  use,  eqns  ( 1 1—3)  and  (11-4)  are  functions  only  of  the  lattice 
frequency  and  its  derivatives.  Lattice  spacing,  per  se,  has  been  removed  completely  as  an 
independent  variable! 

We  may  compare  the  results  of  eqns  (11-3)  and  (1 1-4)  to  the  triple-summed  force  and 
energy  calculations  given  by  eqns  (4-5)  and  (4-2),  respectively.  We  do  so  in  Figure  10  for  two 
cases  of  a.  In  compression  (Figures  10b  and  lOd),  the  analytical  result  for  F  is  within  2%  of  the 
summed  result  to  a  relative  lattice  spacing  of  0.584  (V7V0  =  0.20)  in  the  case  of  a  =  2,  and  0.578 
(V7V0  =  0.19)  in  the  case  of  a  =  5.  In  expansion  (Figures  10a  and  10c),  the  analytical  forms  do 
an  excellent  job  out  to  X  =  Xslah  and  rapidly  diverge  beyond  that  point.  As  mentioned  earlier,  the 
analytical  forms  are  only  defined  out  to  the  lattice  spacing  where  to  goes  to  zero,  which  is  usually 
rather  close  in  value  to  Xslab.  From  this,  we  may  surmise  that  the  analytical  form  does  an 
excellent  job  of  representing  the  lattice  composed  of  discrete  atoms  by  way  of  a  continuous 
function,  as  long  as  the  granularity  of  the  lattice  (spacing  normalized  by  some  function  of  a)  is 
not  too  large.  The  parameter  Xsmb  is  an  approximate  measure  of  this  granularity,  defining  a  cutoff 
point  beyond  which  the  analytical  result  will  eventually  diverge  from  the  discretized  reality  of 
an  actual  lattice.  To  show  that  the  analytical  equations  above  are  capable  of  capturing  the  highly 
expanded  behavior  of  the  lattice,  when  the  granularity  is  low,  we  provide  Figure  11,  depicting 
the  potential  and  force  for  an  a  =  1.1  lattice.  In  this  case,  the  value  of  Xslab  /Xq  does  not  occur 
until  a  value  in  excess  of  7.3.  Thus,  for  lattice  spacings  smaller  than  this  amount,  as  depicted 
in  the  figure,  the  analytical  forms  (dashed  lines)  do  a  good  job  of  matching  the  summed 
equations  (solid  lines). 

Even  though  the  fits  entailed  in  eqns  (10-1  >-(10-4)  were  employed  in  the  development 
of  eqns  (11-3)  and  (11-4),  the  results  portrayed  in  Figures  10  and  11  employ  the  actual  values 
for  to  and  \|r  arising  from  the  triple  summations  and  not  some  fitted  form  for  their  behavior. 
Thus,  eqns  (1 1-3)  and  (1 1-4)  are  seen  not  to  require  an  artificially  imposed  functional  behavior 
for  \|/  or  to  [i.e.,  the  primary  purpose  of  eqn  (10-1)  was  to  assist  in  the  development  of 
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(a) 


m0 

(b) 

Figure  10.  A  comparison  of  the  lattice’s  energy  potential  and  force  (solid  lines)  to  the 
analytical,  frequency-based  result  (dashed  lines):  (a)  FCC,  a  =  2  (expansion);  (b) 
FCC,  a  =  2  (compression);  (c)  FCC,  a  =  5  (expansion);  (d)  FCC,  a  =  5 
(compression).  Shown  on  (a)  and  (c)  are  lattice  spacings  corresponding  to  \tab 
and  co  =  0. 
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Figure  10.  A  comparison  of  the  lattice’s  energy  potential  and  force  (solid  lines)  to  the 
analytical,  frequency-based  result  (dashed  lines):  (a)  FCC,  a  =  2  (expansion);  (b) 
FCC,  a  =  2  (compression);  (c)  FCC,  a  =  5  (expansion);  (d)  FCC,  a  =  5 
(compression).  Shown  on  (a)  and  (c)  are  lattice  spacings  corresponding  to  Xslab 
and  to  =  0. 
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Figure  11.  A  comparison  of  the  summed  energy  and  lattice  force  (solid  lines)  to  the 
analytical  forms  of  eqns  (11-3)  and  (11-4)  (dashed  lines),  for  the  FCC  case  of 
very  low  a  =  1.1. 


49 


eqn  (1 1—2),  and  not  as  an  end  unto  itself).  To  determine  k,  we  may  add  the  individual 
calculations  of  (3  and  eqns  (9-2)  and  (10-2),  which  has  the  effect  of  defining  k  to  match  the 
curvature  of  F  at  the  reference  state.  When  k  is  computed  in  this  manner,  it  may  be  explicitly 
defined  as 


K  =  irlvn  +  1  -  3/2  •  r0 (dy/dV)0 ]/r0  .  d  i-5) 

Computed  values  for  k  (and  for  comparison,  K)  are  given  in  Table  5,  for  a  variety  of  lattice 
conditions.  Note  that,  were  k  to  take  on  the  value  of  K ,  eqn  (11-1)  would  be  recovered  from 
(1 1-4).  Similarly,  as  to  — » 0,  Ec  does  not,  in  general,  approach  Eh,  but  rather  Eh  ( K/k )2.  Since 
Table  5  shows  the  values  for  K  and  k  converging  at  low  a,  however,  we  might  conclude  that, 
in  the  low-anharmonicity  limit,  the  analytical  form  should  match  the  summed  potential  all  the 
way  to  Ec  =  Eh,  which  corresponds  to  X  — » °°.  We,  in  fact,  noted  this  apparent  trend  in  Figure  11. 

One  last  point  to  note  is  that  eqn  (1 1-3)  actually  provides  an  even  better  fit  to  the  triple- 
summed  lattice-force  equations  than  does  eqn  (9-1),  from  which  it  derives.  Presumably,  the 


Table  5.  Values  for  Parameters2  K  and  k 


a 

K  (K) 

SCb 

BCC 

FCC 

i.i 

1.23  (1.26) 

1.26  (1.27) 

1.26  (1.28) 

1.5 

1.22  (1.30) 

1.30  (1.34) 

1.31  (1.35) 

2 

1.19  (1.36) 

1.29  (1.42) 

1.31  (1.43) 

2.5 

1.16  (1.42) 

1.27  (1.48) 

1.29  (1.49) 

3 

1.13  (1.46) 

1.25  (1.53) 

1.27  (1.53) 

4 

1.09  (1.54) 

1.20  (1.59) 

1.22  (1.59) 

5 

1.05  (1.59) 

1.18  (1.63) 

1.18  (1.63) 

a  Computed  from  principal  stiffnesses  only. 

b  For  simple  cubic  lattice,  computed  from  longitudinal  stiffness  only,  as  SC  lattice  not  stable  in  shear,  along  the  principal 
lattice  directions. 
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approximations  in  the  functional  behavior  of  used  to  derive  eqn(ll-3),  compensated 
somewhat  for  the  deviation  noted  in  eqn  (9-1)  at  the  higher  interaction  decay-rates. 

Since  the  current  modeling  has  T  =  r(V ),  which  results  from  the  fact  that  the  vibrational 
stiffness  was  not  modeled  to  be  temperature  dependent,  the  current  results  may  be  cast  in  the 
Griineisen  equation-of-state  form: 


p\\r  -  E  =  pc\\r  -  Ec  .  (11-6) 

Substituting  eqns(ll-3)  and  (11-4)  into  eqn  (11-6)  gives  the  following  equation  of  state: 

P ¥  -  E  =  B0V0/(r0K)2  {[(co/g)0)k  -  1]  +  k(k-  1)  (co/(o0)K  ln(co/co0)>  .  (11-7) 

12.  Comparison  to  Data 


As  the  first  step  in  testing  the  merits  of  eqn  (11-7),  we  will  choose  to  compute  various 
thermodynamic  curves  from  eqn  (1 1-7)  [or  eqn  (1 1-3)  from  which  it  derives]  and  compare  these 
results  to  experimental  data.  Since  these  equations  are  not  a  function  of  X  or  V,  we  must, 
somehow,  relate  them  back  to  lattice  spacing  in  order  to  actually  accomplish  the  comparison. 
The  simplest  approach  would  be  to  assume  the  validity  of  eqn  (10-3),  which  amounts  to 
assuming  a  functional  form  for  T.  The  real  behavior  of  \|r  is,  no  doubt,  more  complicated  than 
this  polynomial  form,  and  one  could  expect  to  achieve  improved  comparisons  to  data  by 
customizing  a  fit  to  for  each  given  case.  Nonetheless,  rather  than  adopting  a  customized  \)r 
function,  the  assumption  of  eqn  (10-3),  along  with  a  fit  to  the  single  parameter  should  provide 
immediate  feedback  on  whether  the  modeling  approach  is  on  the  right  track. 

For  the  case  already  considered,  that  of  specifying  a  pairwise  potential  for  a  lattice,  all 
of  the  parameters  required  to  evaluate  the  merits  of  eqn  (11-7)  are  available  from  the  triple¬ 
summation  equations  for  the  lattice  energy,  force,  and  stiffnesses.  To  evaluate  eqn  (1 1-7)  for  an 
actual  crystalline  material  where  the  pairwise  potential  is  not  known,  however,  requires  the 
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acquisition  of  the  various  model  parameters  by  other  means.  Required  of  the  model  are  values 
for  a  number  of  parameters  at  the  equilibrium  state.  Many  of  them,  such  as  B0  and  V0,  are 
handbook  values  [furthermore,  only  their  product,  defining  the  bulk  sound  speed  via 
C0  =  (B0  V0)1/2,  appears  in  the  model].  Others,  such  as  ro,  are  widely  reported  in  the  literature  as 
well,  though  these  reported  values  may  vary  somewhat  from  source  to  source.  For  the  volumetric 
anharmonicity,  rm/0,  we  use  the  approximate  rule  of  thumb  (reported  here  in  section  8)  that 
rvo/0  =  T|/2.6,  since  values  for  T|  are  provided  for  many  materials  by  Rose  etal.[  11].  In  the  future, 
more  precise  methods  may  be  employed  to  directly  measure  the  curvature  of  the  equilibrium 
lattice  force.  Finally,  in  order  to  obtain  k  by  way  of  eqn  (1 1-5),  a  value  of  (d\\r/dV}0  is  required. 
Though  such  a  value  is  obtainable  from  material  conditions  at  the  equilibrium  state,  there  has 
been  little  or  no  occasion  to  report  such  data  in  the  literature.  Therefore,  we  choose  to  estimate 
this  parameter  by  selecting  a  value  of  ^  (hopefully)  in  the  vicinity  of  0.5,  and  employing 
eqn  (10-2)  to  evaluate  (dy/dV)0  as 


(dy\r/dV)0  =  4/(3  T0)  -  £,  .  (12-1) 

With  all  the  parameters  required  of  the  model  now  obtainable,  we  may  proceed  with  a  few 
representative  comparisons.  We  will  accomplish  the  comparisons  by  calculating  both  the  cold- 
compression  curve  and  shock  Hugoniot  for  several  materials.  The  cold-compression  curve  is  a 
direct  test  of  eqn  (1 1-3)  in  the  absence  of  thermal  effects.  The  shock  Hugoniot,  by  contrast,  is 
a  function  of  both  the  mechanical  and  the  thermal  properties  of  the  lattice.  It  may  be  obtained 
by  eliminating  E  from  eqn  (11-7)  by  way  of  the  Rankine-Hugoniot  shock-energy  relation 

E  -  E0  =  {p0  +  p)(V0  -  V)  /  2  ,  (12-2) 

which  governs  internal  energy  across  a  shock  wave.  For  the  reference  Hugoniot  originating  from 
the  equilibrium  lattice  condition,  both  E0  and  p0  are  zero. 

The  quality  of  this  match  will  test  not  only  eqn  (9-1)  but  also  the  observed  correlation 
of  eqn (10-1),  both  of  which  have  fed  into  the  formulation  of  eqns(ll-3),  (11-4),  and. 


52 


ultimately,  (11-7).  Comparisons  to  silver,  aluminum,  copper,  and  stainless  steel  are  shown  in 
Figures  12-15,  with  the  corresponding  parameters  given  in  Table  6.  The  comparisons  to  both 
the  cold-compression  and  shock-Hugoniot  data  for  these  representative  materials  are  very  good 
to  a  number  of  megabars  of  compression  before,  in  some  cases,  diverging.  And  as  said 
previously,  the  application  of  a  more  general  form  on  \| /,  in  preference  to  the  assumed  form  of 
eqn  (10-3),  can  further  improve  the  correlation  to  data.  The  purpose,  however,  in  retaining  the 
simplified  assumption  of  eqn  (10-3)  is  merely  to  demonstrate  the  general  validity  of  the  method 
and  approach,  as  laid  out  in  this  paper.  Further  refinements  to  the  model,  with  the  intent  to 
improve  correlations  to  data,  would,  at  this  time,  serve  only  to  obfuscate  the  simplicity  of  this 
introductory  report,  which  lays  out  a  static-atomic  paradigm  for  lattice  vibration.  Perhaps  of 
more  importance,  however,  is  the  fact  that  above  several  (1  to  2)  megabars  of  pressure,  Mitchell 
et  al.  [26]  note  that  shocked  metals  are  prone  to  melting.  We  might  therefore  expect  that  a 
model,  such  as  the  present  one,  which  presumes  the  material  to  remain  in  the  solid  state,  might 
no  longer  calibrate  to  the  very  high  shock  pressure  data. 


Table  6.  Parameters  for  Experimental  Comparison 


C0  =  (B0  V0)m 

(m/s) 

1/Vo 

(kg/m3) 

ro 

r  a 

1  voK) 

i 

Ag 

3221 

10490 

2.22 

2.29 

0.500 

1.33 

0.10 

Al 

5189 

2700 

2.03 

1.84 

0.716 

1.49 

-0.06 

Cu 

3995 

8930 

2.02 

2.10 

0.520 

St.  Steel 

4571 

7896 

1.81 

2.00 

0.550 

1.38 

0.19 

a  Selected  in  vicinity  of  =  t|/2.6.  Values  for  T|  obtained  from  Rose  et  al.  [11]. 
b  Computed  by  way  of  eqn  (1 1-5). 
c  Computed  by  way  of  eqn  (12-1). 


We  take  brief  pause,  however,  to  consider  one  difference  noted  between  the  analysis  and 
the  experimental  data — and  that  is  the  generally  low  values  of  T0  computed  in  Table  3.  An 
examination  of  the  modes  constituting  r,  as  shown  in  Figure  5,  can  shed  some  light  on  this 
subject.  From  this  figure,  we  note  that  the  transverse  mode  is  more  stable  than  the  longitudinal 
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v/v0 

(a) 


v/v0 

(b) 

Figure  12.  The  cold-compression  and  shock-Hugoniot  curves  for  silver:  (a)  to 
1.8  megabars;  (b)  to  5  megabars.  Note  that  cold-compression  data  [21]  are 
filled  symbols  and  Hugoniot  data  [22, 23]  are  open  symbols. 
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.5  .6  .7  .8  .9  1.0 


v/v0 

(a) 


.50  .75  1.00 


v/v0 

(b) 

Figure  13.  The  cold-compression  and  shock-Hugoniot  curves  for  aluminum:  (a)  to 
2.5  megabars;  (b)  to  11  megabars.  Note  that  cold-compression  data  [16,24,25] 
are  filled  symbols  and  Hugoniot  data  [23, 26]  are  open  symbols. 
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v/v0 

(b) 

Figure  14.  The  cold-compression  and  shock-Hugoniot  curves  for  copper:  (a)  to 
8  megabars;  (b)  to  25  megabars.  Note  that  cold-compression  data  [21,24]  are 
filled  symbols  and  Hugoniot  data  [23, 26]  are  open  symbols. 
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v/v0 


Figure  15.  The  cold-compression-  and  shock-Hugoniot  curves  for  stainless  steel  to 
4  megabars.  Note  that  cold-compression  data  [11]  are  filled  symbols  and 
Hugoniot  data  [23]  are  open  symbols. 
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one,  in  that  \}/y  remains  defined  to  a  larger  X  than  does  \\rL.  Let  us  keep  in  mind,  here,  that  the 
L  and  T  modes  employed  in  this  analysis  represent  but  three  of  the  many  frequency  modes  of  the 
lattice  and  that,  in  reality,  the  true  \jr  will  be  composed  of  all  the  frequency  modes  and  not  just 
these  three. 

The  makeup  of  the  aggregation  relation,  eqn  (7-1),  is  such  that  as  soon  as  any  vibrational 
mode  goes  unstable  (O),  — >  0),  the  aggregated  (0  must  do  so  as  well.  And  though  the  behavior  of 
the  aggregated  \| ;  tracks  the  more  stable  (larger)  \\rT  at  larger  compressions,  we  do,  in  fact,  see 
that  y  tracks  the  less  stable  \\rL  at  expansions  beyond  Xslab,  the  maximum  in  the  \| r  curve.  If  there 
exists  in  the  lattice  a  mode  less  stable  than  the  principal  L  mode  considered  in  this  report,  then 
the  aggregated  \\r  should  instead  track  that  least-stable  mode  as  the  lattice  expansion  approaches 
the  stability  limit  of  to— >0.  The  net  effect  of  the  modal  coupling  of  eqn (7-1)  is  that,  when 
moving  from  a  compressed  lattice  state  to  one  of  expansion,  begins  to  deviate  from  the 
component  \| r  associated  with  the  most-stable  mode  of  the  lattice,  and  asymptotes  to  the  behavior 
of  the  least-stable  lattice  mode. 

The  lowering  of  \|r  with  increased  X,  arising  from  the  influence  of  the  least  stable  lattice 
mode,  amounts  to  a  net  increase  in  T  over  what  the  stable  mode  would  have  established  (since 
r  =  VA y).  The  fact  that  values  for  T0,  as  calculated  in  Table  3,  are  lower  than  comparable 
experimental  values  allows  us  to  conclude  that  the  principal  L  vibrational  mode  considered  in  this 
study  was  not  the  least- stable  mode  of  the  lattice.  Whereas  we  see  little  deviation  between  \\r  and 
i|/y,  in  Figure  5,  at  the  X0  reference  condition  (\| rL  only  exerts  profound  influence  at  larger 
expansions),  the  higher  values  of  T  provided  by  experimental  data  suggest  that  the  least-stable 
mode,  in  reality,  is  already  influencing  y  at  the  reference  state,  X^.  This  brief  analysis  suggests 
that  future  efforts  along  this  line  of  attack  might  provide  improved  results  by  examining  the  least 
and  most  stable  lattice  modes,  as  opposed  to  the  modes  associated  with  the  principal  directions 
of  the  lattice  structure. 
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13.  Conclusions 


This  report  details  a  simple,  static-atomic  approach  for  computing  the  characteristic 
vibrational  behavior  of  an  atomic  lattice,  when  the  properties  of  the  pairwise  interaction  are 
given.  At  the  core  of  the  method  is  the  understanding  that,  when  atomic  interactions  are 
accounted  for  beyond  the  extent  of  nearest  neighbors,  the  vibrational  and  volumetric  stiffnesses 
are,  by  necessity,  functionally  distinct.  Such  a  distinction  is  at  odds  with  virtually  all  prior 
analytical  treatments  of  the  subject  [2, 5, 6-9],  perhaps  because  accounting  for  the  interactions  of 
nonnearest  neighbors  in  three  dimensions  has,  before  now,  never  been  achievable  in  a 
straighforward,  simple  manner.  Admittedly,  the  stiffness  is  computed  along  only  the  principal 
axes  of  the  lattice  and  so  anisotropy  in  the  stiffness  is  not  accounted  for.  Nonetheless,  these 
limitations  are  the  same  employed  in  a  number  of  earlier  analyses  on  the  subject  [2, 5] — analyses 
that  have  proved  valuable. 

The  results  of  employing  the  method  are  equations  involving  triple  summations,  which 
are  able  to  provide  the  lattice  potential,  interatomic  force,  volumetric-  and  characteristic- 
vibrational  stiffnesses,  as  well  as  the  Griineisen  function  for  the  lattice — in  essence,  the  complete 
equation  of  state  for  the  material.  An  analysis  was  performed  to  show  how  the  mechanical 
behavior  of  the  lattice  compares  to  the  behavior  of  the  constituent  pairwise  interaction.  Results 
also  confirmed  the  distinct  functional  nature  of  the  volumetric  and  vibrational  stiffnesses  of  the 
lattice. 


By  examining  the  current  results  in  light  of  a  recently  published  semi-empirical  model, 
the  triple-summation  equation  for  lattice  force  could  be  cast  in  a  closed-form  analytical  equation, 
in  terms  of  the  characteristic  vibrational  behavior  of  the  lattice,  and  without  the  use  of  fitting 
parameters.  The  thermal  properties  of  the  computed  lattice  were  shown,  and  an  empirical 
correlation  was  observed  between  the  characteristic  vibrational  frequency  of  the  lattice  and  its 
derivative.  This  correlation  seemed  valid,  over  a  wide  range  of  stiffnesses,  for  all  but  the  largest 
expansions  and  vanishingly  small  compressions.  The  correlation,  when  adopted,  translates  into 
analytically  specifiable  functional  forms  for  T,  \jr,  and  co. 
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When  the  analytical  functional  form  put  forth  for  the  mechanical  force  in  the  lattice  is 
combined  with  the  correlation  found  to  govern  the  thermal  properties  of  the  lattice,  a 
comprehensive,  thermal  equation  of  state  is  established.  This  equation,  it  turns  out,  is  identical 
in  form  to  the  recently  proposed  semi-empirical  form  [10],  but  for  the  value  of  a  single  exponent. 
Because  the  form  of  the  equation  of  state  is  in  terms  of  the  characteristic  vibrational  frequency 
of  the  lattice,  this  analytical  equation  of  state  is  only  valid  at  states  where  the  characteristic 
frequency  is  defined.  The  current  analysis  has  determined  that,  at  some  level  of  expansion 
beyond  the  inflection  point  of  the  lattice  potential,  the  vibrational  stiffness  changes  sign  and,  thus, 
the  characteristic  vibrational  frequency  becomes  undefined.  Therefore,  the  analytical  equation 
of  state  proposed  will  be  unable  to  make  thermodynamic  predictions  for  lattice  separations 
beyond  this  critical  value  of  lattice  spacing. 

A  comparison  of  the  analytical  equation  of  state  was  performed  against  cold-compression 
and  shock-Hugoniot  data  for  several  materials.  Without  a  given  pairwise  interaction  known  to 
govern  the  behavior  of  the  lattice,  and  because  the  model  is  a  function  only  of  the  lattice’s 
vibrational  frequency  and  its  derivatives,  an  assumption  had  to  be  made  about  the  the  manner  in 
which  the  vibrational  frequency  varies  with  lattice  spacing.  The  assumption  adopted  was  the 
simplest  possible — namely,  that  the  vibrational  behavior  of  the  actual  lattice  was  identical  in  form 
to  that  noted  (in  the  pairwise-interaction  study)  by  the  correlation  of  eqn  ( 10-1 ).  That  assumption 
requires  the  specification  of  a  single  scalar  parameter,  which  is  known  to  lie  in  the  vicinity  of 
1/2  (the  low-anharmonicity  limit  for  the  parameter).  Though  better  correlations  with  data  are 
obtainable  by  deviating  from  this  assumption,  the  assumption  was  nonetheless  retained,  to 
highlight  both  the  simplicity  of  the  model  and  its  ability  to  predict  the  compressive  and  thermal 
behavior  of  an  actual  lattice  over  a  wide  change  in  volume,  corresponding  to  several  megabars 
of  pressure. 
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NO.  OF 

COPIES  ORGANIZATION 

4  COMMANDER 

US  ARMY  BELVOIR  RD&E  CTR 
ATTN  STRBE  NAE  B  WESTLICH 
STRBE  JMC  T  HANSHAW 
STRBE  NAN 
S  G  BISHOP 
J  WILLIAMS 

FORT  BELVOIR  VA  22060-5166 

3  COMMANDER 

US  ARMY  RESEARCH  OFFICE 
ATTN  J  CHANDRA 
K  IYER 
J  BAILEY 
PO  BOX  12211 

RESEARCH  TRIANGLE  PARK  NC  27709-2211 

1  DIRECTOR 

NAVAL  CIVIL  ENGRNG  LAB 
ATTN  J  YOUNG  CODE  L56 
PORT  HUENEME  CA  93043 

1  NAVAL  POSTGRADUATE  SCHOOL 
PHYSICS  DEPARTMENT 
ATTN  JOSEPH  STERNBERG 
MONTEREY  CA  93943 

1  NAVAL  AIR  WARFARE  CTR 
ATTN  STEPHEN  A  FINNEGAN 
BOX  1018 

RIDGECREST  CA  93556 

3  COMMANDER 

NAVAL  WEAPONS  CENTER 
ATTN  T  T  YEE  CODE  3263 
D  THOMPSON  CODE  3268 
W  J  MCCARTER  CODE  6214 
CHINA  LAKE  CA  93555 
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12  COMMANDER 

NAVAL  SURFACE  WARFARE  CTR 
DAHLGREN  DIVISION 
ATTN  H  CHEN 
D  L  DICKINSON  CODE  G24 
CHARLES  R  ELLINGTON 
C  R  GARRETT  CODE  G22 
W  HOLT  CODE  G22 
R  MCKEOWN 
W  WALLACE  MORTON  JR 
JOHN  M  NELSON 
M  J  SILL  CODE  Hll 
WILLIAM  J  STROTHER 
A  B  WARDLAW,  JR. 

L  F  WILLIAMS  CODE  G33 
17320  DAHLGREN  RD 
DAHLGREN  VA  22448 

5  AIR  FORCE  ARMAMENT  LAB 
ATTN  AFATL  DLJW 
W  COOK 
M  NIXON 

AFATL  DUR  J  FOSTER 
AFATL  MNW 
LT  D  LOREY 
R  D  GUBA 
EGLIN  AFB  FL  32542 

1  USAF  PHILLIPS  LABORATORY 
VTSI 

ATTN  ROBERT  ROYBAL 
KIRTLAND  AFB  NM  87117-7345 

2  USAF  PHILLIPS  LABORATORY 
ATTN  PL  WSCD  FIROOZ  ALLAHDADI 
PV  VTA  DAVID  SPENCER 

3550  ABERDEEN  AVE  SE 
KIRTLAND  AFB  NM  87117-5776 

5  WRIGHT  LABS 

ATTN  MNMW  JOEL  W  HOUSE 

ARMAMENT  DIRECTORATE  STE  326  B1 

RONALD  D  HUNT 

BRYAN  MILLIGAN 

BRUCE  C  PATTERSON 

WADE  H  VAUGHT 

101  W  EGLIN  BLVD 

EGLIN  AFB  FL  32542-6810 


1  AFIT  ENC 

ATTN  DAVID  A  FULK 

WRIGHT  PATTERSON  AFB  OH  45433 

42  DIRECTOR 
LANL 

ATTN  M  LUCERO  MS  A105 
J  V  REPA  MS  A133 
J  P  RITCHIE  MS  B214  T14 
R  DINGUS  MS  B218 
N  KRKORIAN  MS  B228 
R  KIRKPATRICK  MS  B229 
R  THURSTON  MS  B229 
C  T  KLINGNER  MS  B294 
R  MILLER  MS  B294 
B  SHAFER  MS  C931 
G  GISLER  MS  D436 
C  RAGAN  MS  D449 
B  LAUBSCHER  MS  D460 
M  O  SCHNICK  MS  F607 
R  WELLS  MS  F607 
R  KOPP  MS  F645 
R  STELLINGWERFMS  F645 
C  WINGATE  MS  F645 
T  ADAMS  MS  F663 
R  GODWIN  MS  F663 
K  JACOBY  MS  F663 
W  SPARKS  MS  F663 
E  J  CHAPYAK  MS  F664 
J  SHANER  MS  F670 
G  CANA  VAN  MS  F675 
R  GREINER  MS  G740 
J  HILLS  MS  G770 
B  HOGAN  MS  G770 
J  BOLSTAD  MS  G787 
J  WALSH  MS  G787 
R  DAVIDSON  MS  K557 
R  HENNINGER  MS  K557  N6 
T  ROLLET  MS  K574 
PHOWE  MSP915 
W  DEAL  MS  P915 
J  KENNEDY  MS  P915 
A  ROACH  MS  P915 
W  HEMSING  MS  P940 
E  POGUE  MS  P940 
J  MCAFEE  MS  P950 
D  PAISLEY  MS  P950 
L  PICKLESIMER  MS  P950 
PO  BOX  1663 
LOS  ALAMOS  NM  87545 
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6  DIRECTOR 
LANL 

ATTN  R  WARNES  MS  P950 
S  SHEFFIELD  MS  P952 
D  MANDELL 
K  MARK 
S  J  MOSSO 
L  SCHWALBE 
PO  BOX  1663 
LOS  ALAMOS  NM  87545 

35  DIRECTOR 

SANDIA  NATL  LABS 
ATTN  E  H  BARSISMS-03 1 
ERIC  W  REECE  MS-0307 
DANIEL  P  KELLYMS-0307 
L  WEIRICK  MS-0327 
R  TACHAU  MS-0425 
D  LONGCOPE  MS-0439 
D  HAYES  MS-0457 
J  ASAY  MS-0458 
W  TEDESCHI  MS-0482 
J  SCHULZE  MS-0483 
P  A  LONGMIREMS-0560 
J  COREY  MS-0576 
E  S  HERTEL  JRMS-0819 
A  ROBINSON  MS-0819 
T  TRUCANO  MS-0819 
J  M  MCGLAUN  MS-0819 
R  BRANNON  MS-0820 
L  CHHABILDASMS-0821 
J  ANG  MS-0821 
M  BOSLOUGH  MS-0821 
L  CHHABILDAS  MS-0821 
D  CRAWFORD  MS-0821 
M  FURNISH  MS-0821 
CHALL  MS-0821 
W  REINHART  MS-0821 
P  STANTON  MS-0821 
M  KIPP  DIV  1533 
P  YARRINGTON  DIV  1533 
J  MCGLAWA  DIV  1541 
M  FORRESTAL  DIV  1551 
R  LAFARGE  DIV  1551 
C  HILLS  DIV  1822 
R  O  NELLUMS  DIV  9122 
P  TAYLOR  ORG  1432 
D  KERNAN  ORG  1433 
PO  BOX  5800 

ALBUQUERQUE  NM  87185 


6  DIRECTOR 
SANDIA  NATL  LABS 
ATTN  B  LEVIN  ORG  7816 
L  N  KMETYK 

R  REEDER 
J  SOUTHWARD 
C  KONRAD 
K  LANG 
PO  BOX  5800 

ALBUQUERQUE  NM  87185 

3  DIRECTOR 
LLNL 
MS  L35 

ATTN  R  E  TIPTON 
D  BAUM 
T  MCABEE 
PO  BOX  808 
LIVERMORE  CA  94550 

7  DIRECTOR 
LLNL 

MS  L122 
ATTN  R  PIERCE 
R  ROSINKY 
O  J  ALFORD 
D  STEWART 
T  VIDLAK 
B  R  BOWMAN 
W  DIXON 
PO  BOX  808 
LIVERMORE  CA  94550 

2  DIRECTOR 
LLNL 
MS  L125 

ATTN  DOUGLAS  R  FAUX 
NORMAN  W  KLINO 
PO  BOX  808 
LIVERMORE  CA  94550 

1  DIRECTOR 
LLNL 

ATTN  ROBERT  BARKER  LI 59 
PO  BOX  808 
LIVERMORE  CA  94550 
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3  DIRECTOR 
LLNL 
MS  L163 

ATTN  MILTON  FINGER 
R  PERRET 
W  SHOTTS 
PO  BOX  808 
LIVERMORE  CA  94550 

3  DIRECTOR 
LLNL 
MS  L178 

ATTN  H  KRUGER 
G  POMYKAL 

MICHAEL  GERASSIMENKO 
PO  BOX  808 
LIVERMORE  CA  94550 

2  DIRECTOR 
LLNL 
MS  L180 

ATTN  G  SIMONSON 
A  SPERO 
PO  BOX  808 
LIVERMORE  CA  94550 

1  DIRECTOR 
LLNL 

ATTN  FRANK  A  HANDLER  LI  82 
PO  BOX  808 
LIVERMORE  CA  94550 

2  DIRECTOR 
LLNL 

MS  L282 
ATTN  W  TAO 
P  URTIEW 
PO  BOX  808 
LIVERMORE  CA  94550 

2  DIRECTOR 
LLNL 
MS  L290 
ATTN  A  HOLT 
J  E  REAUGH 
PO  BOX  808 
LIVERMORE  CA  94550 


1  DIRECTOR 
LLNL 

ATTN  W  J  NELLIS  L299 
PO  BOX  808 
LIVERMORE  CA  94550 

1  DIRECTOR 
LLNL 

ATTN  D  WOOD  L352 
PO  BOX  808 
LIVERMORE  CA  94550 

1  DIRECTOR 
LLNL 

ATTN  STEPHEN  G  COCHRAN  L389 
PO  BOX  808 
LIVERMORE  CA  94550 

3  DIRECTOR 
LLNL 
MS  L495 
ATTN  D  GAVEL 
J  HUNTER 
E  JOHANSSON 
PO  BOX  808 
LIVERMORE  CA  94550 

1  DIRECTOR 
LLNL 

ATTN  R  M  KUKLO  L874 
PO  BOX  808 
LIVERMORE  CA  94550 

8  DIRECTOR 
LLNL 

ATTN  G  W  REPP 
M  J  MURPHY 
A  C  MITCHELL 
J  A  MORI  ARTY 
R  A  HEINLE 
N  C  HOLMES 
M  SHANNON 
BMDO  ROBERT  M  HALL 
PO  BOX  808 
LIVERMORE  CA  94550 
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1 
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ENERGETIC  MATERIALS  RSCH  CTR/DOE 

NEW  MEXICO  INST  OF  MINING  &  TECH 

ATTN  DAVID  J  CHAVEZ 

LARRY  LIBERSKY 

FRED  SANDSTROM 

CAMPUS  STATION 

SOCORRO  NM  87801 

NASA 

LEWIS  RESEARCH  CENTER 
ATTN  J  FERRANTE 
CLEVELAND  OH  44135 

NASA 

JOHNSON  SPACE  CENTER 
ATTN  ERIC  CHRISTIANSEN 
JEANNE  LEE  CREWS 
FREDRICH  HORZ 
MAIL  CODE  SN3 
2101  NASA  RD  I 
HOUSTON  TX  77058 

APPLIED  RESEARCH  LAB 
ATTN  JEFFREY  A  COOK 
10000  BURNETT  ROAD 
AUSTIN  TX  78758 

GM  RESEARCH  LABS 
ATTN  J  R  SMITH 
WARREN  MI  48090 

JET  PROPULSION  LABORATORY 
IMPACT  PHYSICS  GROUP 
ATTN  ZDENEK  SEKANINA 
PAUL  WEISSMAN 
BOB  WEST 
JAMES  ZWISSLER 
MARC  ADAMS 
4800  OAK  GROVE  DR 
PASADENA  CA  91 109 

MIT  LINCOLN  LAB 
ARMY  SCIENCE  BOARD 
ATTN  WADE  M  KORNEGAY 
244  WOOD  ST  RM  S2  139 
LEXINGTON  MA  02173 


1  BOSTON  UNIVERSITY 
DEPT  OF  PHYSICS 
ATTN  ZEEV  JAEGER 

590  COMMONWEALTH  AVE 
BOSTON,  MA  02215 

2  BROWN  UNIVERSITY 
ATTN  R  CLIFTON  (ENGNG) 

P  SCHULTZ  (GEO  SCI) 

PROVIDENCE  RI  02912 

3  CALTECH 

ATTN  J  SHEPHERD  MS  105-50 
ANDREW  P  INGERSOLL  MS  170-25 
THOMAS  J  AHRENS  MS  252-21 
1201  E  CALIFORNIA  BLVD 
PASADENA  CA  91125 

1  CALTECH 

ATTN  GLENN  ORTON  MS  169  237 
4800  OAK  GROVE  DR 
PASADENA  CA  91007 

3  CORNELL  UNIVERSITY 

DEPT.  MATERIALS  SCIENCE  &  ENGNG 
ATTN  R  G  GREENE 
H  LUO 
A  L  RUOFF 
ITHACA  NY  14853 

3  DREXEL  UNIVERSITY 
ATTN  MEM  DEPT 
PHYSICS  DEPT 

A  ZAVALIANGOS  (DEPT  MAT  ENGNG) 
32ND  &  CHESTNUT  ST 
PHILADELPHIA  PA  19104 

1  GEORGIA  INSTITUTE  OF  TECHNOLOGY 
COMPUTATIONAL  MODELING  CENTER 
ATTN  S  ATLURI 
ATLANTA  GA  30332-0356 

1  GEORGIA  INSTITUTE  OF  TECHNOLOGY 
SCHOOL  OF  MATL  SCIENCE  &  ENGNG 
ATTN  K  LOGAN 
ATLANTA  GA  30332-0245 
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1 

IOWA  STATE  UNIVERSITY 

DEPT  PHYSICS  AND  ASTRONOMY 

ATTN  JIM  ROSE 

34  PHYSICS 

AMES  I A  50011 

1 

UC  BERKELEY 

MECHANICAL  ENGINEERING  DEPT 
GRADUATE  OFFICE 

ATTN  KEZHUN  LI 

BERKELEY  CA  94720 

1 

JOHNS  HOPKINS  UNIVERSITY 

MAT  SCI  &  ENGNG  DEPT 

ATTN:  MO  LI 

102  MARYLAND  HALL 

3400  N  CHARLES  ST 

1 

UC  DAVIS 

INST  OF  THEORETICAL  DYNAMICS 
ATTN  E  G  PUCKETT 

DAVIS  CA  95616 

BALTIMORE  MD  21218-2689 

1 

UC  LOS  ANGELES 

DEPT  OF  MAT  SCIENCE  &  ENGNG 

5 

JOHNS  HOPKINS  UNIVERSITY 

APPLIED  PHYSICS  LAB 

ATTN  TERRY  R  BETZER 

ATTN  J  J  GILMAN 

LOS  ANGELES  CA  90024 

ALVIN  R  EATON 

RICHARD  H  KEITH 

DALE  K  PACE 

ROGER  L  WEST 

JOHNS  HOPKINS  ROAD 

LAUREL  MD  20723 

2 

UC  SAN  DIEGO 

DEPT  APPL  NECH  &  ENGR 

SVCS  R011 

ATTN  S  NEMAT-NASSER 

M  MEYERS 

LA  JOLLA  CA  92093-041 1 

1 

LOUISIANA  STATE  UNIVERSITY 

ATTN  ROBERT  W  COURTER 

948  WYLIE  DR 

BATON  ROUGE  LA  70808 

2 

UNIV  OF  ALA  HUNTSVILLE 
AEROPHYSICS  RSCH  CTR 

ATTN  GARY  HOUGH 

DAVID  J  LIQUORNIK 

PO  BOX  999 

1 

NC  STATE  UNIVERSITY 

ATTN  YASUYUKI  HORIE 

HUNTSVILLE  AL  35899 

RALEIGH  NC  27695-7908 

1 

UNIV  OF  ALA  HUNTSVILLE 

MECH  ENGRNG  DEPT 

1 

PENNSYLVANIA  STATE  UNIVERSITY 
ATTN  PHYSICS  DEPT 

UNIVERSITY  PARK  PA  16802 

1 

ATTN  W  P  SCHONBERG 

HUNTSVILLE  AL  35899 

UNIVERSITY  OF  CHICAGO 

4 

SOUTHWEST  RESEARCH  INSTITUTE 

ATTN  C  ANDERSON 

S  A  MULLIN 

J  RIEGEL 

J  WALKER 

DEPT  OF  THE  GEOPHYSICAL  SCIENCES 
ATTN  G  H  MILLER 

5734  S  ELLIS  AVE 

CHICAGO  IL  60637 

PO  DRAWER  28510 

SAN  ANTONIO  TX  78228-0510 

3 

UNIVERSITY  OF  DAYTON  RSCH  INST 
ATTN  N  BRAR 

D  GROVE 

1 

TEXAS  A&M  UNIVERSITY 

PHYSICS  DEPARTMENT 

ATTN  DAN  BRUTON 

COLLEGE  STATION  TX  77843-4242 

A  PIEKUTOWSKI 

300  COLLEGE  PARK 

DAYTON  OH  45469-0182 
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1 


2 


ORGANIZATION 


NO.  OF 

COPIES  ORGANIZATION 


UNIVERSITY  OF  DELAWARE 

DEPT  OF  MECHANICAL  ENGINEERING 

ATTN  PROF  J  GILLESPIE 

DEAN  R  B  PIPES 

PROF  J  VINSON 

PROF  D  WILKINS 

NEWARK  DE  19716 

UNIVERSITY  OF  ILLINOIS 
PHYSICS  BUILDING 
ATTN  A  V  GRANATO 
URBANA,  IL  61801 

UNIVERSITY  OF  MARYLAND 
ATTN  PHYSICS  DEPT  (BLDG  082) 
COLLEGE  PARK  MD  20742 

UNIVERSITY  OF  PENNSYLVANIA 
ATTN  P  A  HEINEY 
DEPT  OF  PHYSICS  &  ASTRONOMY 
209  SOUTH  33RD  ST 
PHILADELPHIA  PA  19104 

UNIVERSITY  OF  PUERTO  RICO 
DEPT  CHEMICAL  ENGINEERING 
ATTN  L  A  ESTEVEZ 
MAYAGUEZ  PR  00681-5000 

UNIVERSITY  OF  TEXAS 
DEPT  OF  MECHANICAL  ENGINEERING 
ATTN  ERIC  P  FAHRENTHOLD 
AUSTIN  TX  78712 

VIRGINIA  POLYTECHNIC  INSTITUTE 
COLLEGE  OF  ENGINEERING 
ATTN  R  BATRA 
BLACKSBURG  VA  24061-0219 

AEROJET 

ATTN  J  CARLEONE 
S  KEY 

PO  BOX  13222 

SACRAMENTO  CA  95813-6000 

AEROJET  ORDNANCE 
ATTN  P  WOLF 
G  PADGETT 
1100  BULLOCH  BLVD 
SOCORRO  NM  87801 


2  ALL1ANT  TECHSYSTEMS  INC 
ATTN  R  STRYK 
G  R  JOHNSON  MNI 1-1614 
600  SECOND  ST  NE 
HOPKINS  MN  55343 

1  MARVIN  L  ALME 

2180  LOMA  LINDA  DR 
LOS  ALAMOS  NM  87544-2769 

1  APPLIED  RESEARCH  ASSOC  INC 
ATTN  JEROME  D  YATTEAU 

5941  S  MIDDLEFIELD  RD  SUITE  100 
LITTLETON  CO  80123 

2  APPLIED  RESEARCH  ASSOC  INC 
ATTN  DENNIS  GRADY 
FRANK  MAESTAS 

SUITE  A220 

4300  SAN  MATEO  BLVD  NE 
ALBUQUERQUE  NM  87110 

1  ATA  ASSOCIATES 
ATTN  W  ISBELL 
PO  BOX  6570 

SANTA  BARBARA  CA  931 1 1 

1  BATTELLE 

ATTN  ROBER  M  DUGAS 

7501  S  MEMORIAL  PKWY  SUITE  101 

HUNTSVILLE  AL  35802-2258 

3  BOEING  AEROSPACE  CO 

SHOCK  PHYSICS  &  APPLIED  MATH 
ENGINEERING  TECHNOLOGY 
ATTN  R  HELZER 
T  MURRAY 
J  SHRADER 
PO  BOX  3999 
SEATTLE  WA  98124 

1  BOEING  HOUSTON  SPACE  STN 
ATTN  RUSSELL  F  GRAVES 
BOX  58747 
HOUSTON  TX  77258 

1  BRIGS  CO 

ATTN  JOSEPH  E  BACKOFEN 
2668  PETERSBOROUGH  ST 
HERNDON  VA  20171-2443 
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1 

CALIFORNIA  RSCH  &  TECHNOLOGY 

ATTN  M  MAJERUS 

PO  BOX  2229 

PRINCETON  NJ  08543 

3 

DYNASEN 

ATTN  JACQUES  CHAREST 

MICHAEL  CHAREST 

MARTIN  LILLY 

20  ARNOLD  PL 

1 

CENTURY  DYNAMICS  INC 

ATTN  N  BIRNBAUM 

GOLETA  C A  93117 

1 

2333  SAN  RAMON  VALLEY  BLVD 

SAN  RAMON  CA  94583-1613 

COMPUTATIONAL  MECHANICS 

1 

R  J  EICHELBERGER 

409  W  CATHERINE  ST 

BEL  AIR  MD  21014-3613 

CONSULTANTS 

ATTN  J  A  ZUKAS 

PO  BOX  11314 

BALTIMORE  MD  21239-0314 

1 

ELORET  INSTITUTE 

ATTN  DAVID  W  BOGDANOFF  MS  230  2 
NASA  AMES  RESEARCH  CENTER 
MOFFETT  HELD  CA  94035 

1 

CYPRESS  INTERNATIONAL 

ATTN  A  CAPONECCHI 

1201  E  ABINGDON  DR 

ALEXANDRIA  VA  22314 

3 

ENIG  ASSOCIATES  INC 

ATTN  J  ENIG 

D  J  PASTINE 

M  CO WPERTHW AITE 

SUITE  500 

1 

DEFENSE  TECHNOLOGY  INTL.  INC 

ATTN  D  E  AYER 

THE  STARK  HOUSE 

11120  NEW  HAMPSHIRE  AVE 

SILVER  SPRING  MD  20904-2633 

22  CONCORD  ST 

NASHUA  NH  03060 

1 

EXPLOSIVE  TECHNOLOGY 

ATTN  M  L  KNAEBEL 

PO  BOX  KK 

1 

DESKIN  RESEARCH  GROUP  INC 

ATTN  EDWARD  COLLINS 

FAIRFIELD  CA  94533 

2270  AGNEW  RD 

SANTA  CLARA  CA  95054 

1 

GB  TECH  LOCKHEED 

ATTN  JAY  LAUGHMAN 

2200  SPACE  PARK  SUITE  400 

3 

DOW  CHEMICAL  INC 

ORDNANCE  SYSTEMS 

HOUSTON  TX  77258 

ATTN  C  HANEY 

A  HART 

B  RAFANIELLO 

800  BUILDING 

MIDLAND  MI  48667 

2 

GB  TECH  LOCKHEED 

ATTN  LUCILLE  BORREGO  C23C 

JOE  FALCON  JR  C23C 

2400  NASA  ROAD  1 

HOUSTON  TX  77058 

1 

G  E  DUVALL 

5814  NE  82ND  COURT 

VANCOUVER  WA  98662-5944 

6 

GDLS 

38500  MOUND  RD 

ATTN  W  BURKE  MZ436-21-24 

G  CAMPBELL  MZ436-30-44 

3 

DYNA  EAST  CORP 

ATTN  P  C  CHOU 

R  CICCARELLI 

W  FLIS 

3620  HORIZON  DRIVE 

KING  OF  PRUSSIA  PA  19406 

D  DEBUSSCHER  MZ436-20-29 

J  ERIDON  MZ436-21-24 

W  HERMAN  MZ  435-01-24 

S  PENTESCU  MZ436-21-24 

STERLING  HTS  MI  48310-3200 
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2  GENERAL  RESEARCH  CORP 
ATTN  A  CHARTERS 
T  MENNA 
PO  BOX  6770 

SANTA  BARBARA  CA  93160-6770 

1  GRC  INTERNATIONAL 

ATTN  TIMOTHY  M  CUNNINGHAM 
5383  HOLLISTER  AVE 
SANTA  BARBARA  CA  931 1 1 

7  INST  OF  ADVANCED  TECHNOLOGY 
UNIVERSITY  OF  TX  AUSTIN 
ATTN  S  J  BLESS 
J  CAZAMIAS 
J  DAVIS 
H  D  FAIR 
T  M  K.IEHNE 
D  LITTLEFIELD 
M  NORMANDIA 
4030-2  W  BRAKER  LN 
AUSTIN  TX  78759 

1  INTERNATIONAL  RESEARCH  ASSOC 
ATTN  D  L  ORPHAL 
4450  BLACK  AVE 
PLEASANTON  CA  94566 

1  INTERPLAY 

ATTN  F  E  WALKER 
18  SHADOW  OAK  RD 
DANVILLE  CA  94526 

1  ITT  SCIENCES  AND  SYSTEMS 
ATTN  J  WILBECK 
600  BLVD  SOUTH,  SUITE  208 
HUNTSVILLE  AL  35802 

1  R  JAMESON 
624  ROWE  DR 
ABERDEEN  MD  21001 

1  KAMAN  SCIENCES  CORP 
ATTN  DENNIS  L  JONES 
2560  HUNTINGTON  AVE  SUITE  200 
ALEXANDRIA  VA  22303 


7  KAMAN  SCIENCES  CORP 
ATTN  J  ELDER 
RICHARD  P  HENDERSON 
DAVID  A  PYLES 
FRANK  R  SAVAGE 
JAMES  A  SUMMERS 
TIMOTHY  W  MOORE 
THY  YEM 

600  BLVD  S  SUITE  208 
HUNTSVILLE  AL  35802 

3  KAMAN  SCIENCES  CORP 
ATTN  SHELDON  JONES 
GARY  L  PADEREWSKI 
ROBERT  G  PONZINI 

1500  GRDN  OF  THE  GODS  RD 
COLORADO  SPRINGS  CO  80907 

4  KAMAN  SCIENCES  CORP 
ATTN  NASIT  ARI 
STEVE  R  DIEHL 
WILLIAM  DOANE 
VERNON  M  SMITH 

PO  BOX  7463 

COLORADO  SPRINGS  CO  80933-7463 

1  D  R  KENNEDY  &  ASSOC  INC 
ATTN  D  KENNEDY 
PO  BOX  4003 

MOUNTAIN  VIEW  CA  94040 

1  KERLEY  PUBLISHING  SERVICES 
ATTN  G  I  KERLEY 

PO  BOX  13835 

ALBUQUERQUE  NM  87192-3835 

2  KTECH  CORPORATION 
ATTN  FRANK  W  DAVIES 
LARRY  M  LEE 

901  PENNSYLVANIA  NE 
ALBUQUERQUE  NM  87110 

1  LIVERMORE  SOFTWARE  TECH  CORP 
ATTN  J  O  HALLQUIST 
2876  WAVERLY  WAY 
LIVERMORE  CA  94550 
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1  LOCKHEED  MARTIN  MISSLE  &  SPACE 
ATTN  WILLIAM  R  EBERLE 
PO  BOX  070017 
HUNTSVILLE  AL  35807 

3  LOCKHEED  MARTIN  MISSILE  &  SPACE 
ATTN  M  A  LEVIN  ORG  81  06  BLDG  598 
M  R  MCHENRY 

T  A  NGO  ORG  81  10  BLDG  157 
111  LOCKHEED  WAY 
SUNNYVALE  CA  94088 

4  LOCKHEED  MISSILE  &  SPACE  CO 
ATTN  JOHN  R  ANDERSON 
WILLIAM  C  KNUDSON 

S  KUSUMI  0  81  11  BLDG  157 
J  PHILLIPS  0  54  50 
PO  BOX  3504 
SUNNYVALE  CA  94088 

1  LOCKHEED  MISSILE  &  SPACE  CO 
ATTN  R  HOFFMAN 
SANTA  CRUZ  FACILITY 
EMPIRE  GRADE  RD 
SANTA  CRUZ  CA  95060 

1  LOCKHEED  NASA  JSC 

SPACE  SCIENCE  BRANCH 
ATTN  JAMES  HYDE 
BOX  58561  MC  B22 
HOUSTON  TX  77258 

1  MCDONNELL  DOUGLAS 
ASTRONAUTICS  CO 
ATTN  B  L  COOPER 
5301  BOLSA  AVE 
HUNTINGTON  BEACH  CA  92647 

2  ORLANDO  TECHNOLOGY  INC 
ATTN  DANIEL  A  MATUSKA 
MICHAEL  GUNGER 

PO  BOX  855 

SHALIMAR  FL  32579-0855 

1  PHYSICAL  SCIENCES  INC 
ATTN  PETER  NEBOLSINE 
20  NEW  ENGLAND  BUS  CTR 
ANDOVER  MA  01810 


3  PHYSICS  INTERNATIONAL 
ATTN  R  FUNSTON 
G  FRAZIER 
L  GARNETT 
PO  BOX  5010 
SAN  LEANDRO  CA  94577 

1  PRC  INC 

ATTN  J  ADAMS 

5166  POTOMAC  DR  #103 

KING  GEORGE  VA  22485-5824 

1  RAYTHEON  ELECTRONIC  SYSTEMS 
ATTN  R  KARPP 
50  APPLE  HILL  DRIVE 
TEWKSBURY  MA  01876 

1  ROCKWELL  INTERNATIONAL 
ROCKETDYNE  DIVISION 
ATTN  H  LEIFER 
16557  PARK  LN  CIRCLE 
LOS  ANGELES  CA  90049 

1  ROCKWELL  MISSILE  SYS  DIV 
ATTN  T  NEUHART 
1800  SATELLITE  BLVD 
DULUTH  GA  30136 

1  SAIC 

ATTN  MICHAEL  W  MCKAY 
10260  CAMPUS  POINT  DR 
SAN  DIEGO  CA  92121 

1  SHOCK  TRANSIENTS  INC 
ATTN  DAVID  DAVISON 
BOX  5357 

HOPKINS  MN  55343 

2  SIMULATION  &  ENG  CO  INC 
ATTN  ELSA  I  MULLINS 
STEVEN  E  MULLINS 

8840  HWY  20  SUITE  200  N 
MADISON  AL  35758 

2  SOUTHERN  RESEARCH  INSTITUTE 
ATTN  LINDSEY  A  DECKARD 
DONALD  P  SEGERS 
PO  BOX  55305 

BIRMINGHAM  AL  35255-5305 
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5  SRI  INTERNATIONAL 
ATTN  JAMES  D  COLTON 
D  CURRAN 
R  KLOOP 
R  L  SEAMAN 
D  A  SHOCKEY 
333  RAVENSWOOD  AVE 
MENLO  PARK  CA  94025 

2  TELEDYNE  BROWN  ENGR 
ATTN  JIM  W  BOOTH 
MARTIN  B  RICHARDSON 
PO  BOX  070007  MS  50 
HUNTSVILLE  AL  35807-7007 

1  ZERNOW  TECHNICAL  SVCS  INC 
ATTN  LOUIS  ZERNOW 
425  W  BONITA  AVE  SUITE  208 
SAN  DIMAS  CA  91773 

1  SUNY  STONEYBROOK 

DEPT  APPL.  MATH  &  ST  AT. 
ATTN  J  GLIMM 
STONEYBROOK  NY  11794 


ABERDEEN  PROVING  GROUND 

9  DIR,  USARL 

ATTN:  AMSRL-WM,  I  MAY 
AMSRL-WM-BC,  A  ZIELINSKI 
AMSRL-WM-BD, 

R  PESCE-RODRIGUEZ 
A  KOTLAR 

AMSRL-WM-BE,  S  HOWARD 
AMSRL-WM-MB,  G  GAZONAS 
AMSRL-WM-MC,  J  M  WELLS 
AMSRL-WM-T, 

W  F  MORRISON 
T  W  WRIGHT 
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43  DIR.  USARL 

ATTN:  AMSRL-WM-TA, 

W  GILLICH 
S  BILYK 
M  BURKINS 
W  BRUCHEY 
J  DEHN 
G  FILBEY 
W  A  GOOCH 
H  W  MEYER 
E  J  RAPACKI 
J  RUNYEON 
AMSRL-WM-TB, 

R  FREY 
P  BAKER 
R  LOTTERO 
J  STARKENBERG 
AMSRL-WM-TC, 

W  S  DE  ROSSET 
T  W  BJERKE 
R  COATES 
F  GRACE 
K  KIMSEY 
M  LAMPSON 
D  SCHEFFLER 
S  SCHRAML 
G  S1LSBY 
B  SORENSEN 
R  SUMMERS 
W  WALTERS 
AMSRL-WM-TD, 

S  CHOU 
A  M  DIETRICH 
J  M  BOTELER 
D  DANDEKAR 
K  FRANK 
M  RAFTENBERG 
A  RAJENDRAN 
M  SCHEIDLER 
S  SCHOENFELD 
S  SEGLETES  (5  CP) 
T  WEERASOORIYA 
AMSRL-WM-WD, 

J  POWELL 
A  PRAKASH 
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1 
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1 
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AERONAUTICAL  &  MARITIME 
RESEARCH  LABORATORY 
ATTN  N  BURMAN 
R  WOODWARD 
S  CIMPOERU 
D  PAUL 
PO  BOX  4331 
MELBOURNE  VIC  3001 
AUSTRALIA 

ABTEILUNG  FUER  PHYSIKALISCHE 
CHEMIE 

MONTANUNIVERSITAET 
ATTN  E  KOENIGSBERGER 
A  8700  LEOBEN 
AUSTRIA 

PRB  S  A 

ATTN  M  VANSNICK 

AVENUE  DE  TERVUEREN  168  BTE  7 

BRUSSELS  B  1150 

BELGIUM 

ROYAL  MILITARY  ACADEMY 
ATTN  E  CELENS 
RENAISSANCE  AVE  30 
B1040  BRUSSELS 
BELGIUM 

BULGARIAN  ACADEMY  OF  SCIENCES 
SPACE  RESEARCH  INSTITUTE 
ATTN  VALENTIN  GOSPODINOV 
1000  SOFIA  PO  BOX  799 
BULGARIA 

CANADIAN  ARSENALS  LTD 
ATTN  P  PELLETIER 
5  MONTEE  DES  ARSENAUX 
VILLIE  DE  GRADEUR  PQ  J5Z2 
CANADA 

DEFENCE  RSCH  ESTAB  SUFFEELD 
ATTN  D  MACKAY 

RALSTON  ALBERTA  TOJ  2NO  RALSTON 
CANADA 


1  DEFENCE  RSCH  ESTAB  SUFFIELD 
ATTN  CHRIS  WEICKERT 
BOX  4000  MEDICINE  HAT 
ALBERTA  TIA  8K6 
CANADA 

1  DEFENCE  RSCH  ESTAB  VALC ARTIER 
ARMAMENTS  DIVISION 
ATTN  R  DELAGRAVE 
2459  PIE  XI  BLVD  N 
PO  BOX  8800 

CORCELETTE  QUEBEC  GOA  1R0 
CANADA 

1  UNIVERSITY  OF  GUELPH 
PHYSICS  DEPT 
ATTN  C  G  GRAY 
GUELPH  ONTARIO 
NIG  2W1 
CANADA 

1  CEA 

ATTN  ROGER  CHERET 
CEDEX  15 

313  33  RUE  DE  LA  FEDERATION 

PARIS  75752 

FRANCE 

1  CEA 

CISI  BRANCH 
ATTN  PATRICK  DAVID 
CENTRE  DE  SACLAY  BP  28 
GIF  SUR  YVETTE  91192 
FRANCE 

1  CEA/CESTA 

ATTN  ALAIN  GEILLE 
BOX  2  LE  BARP  33114 
FRANCE 

6  CENTRE  D’ETUDES  DE  GRAMAT 
ATTN  SOLVE  GERARD 
CHRISTIAN  LOUPIAS 
PASCALE  OUTREBON 
J  CAGNOUX 
C  GALLIC 
J  TRANCHET 
GRAMAT  46500 
FRANCE 
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2 

CENTRE  D’ETUDES  DE  LIMEIL- V ALENTON 
ATTN  CHRISTIAN  AUSSOURD 

JEAN-CLAUDE  BOZIER 

SAINT  GEORGES  CEDEX 

VILLENEUVE  94195 

FRANCE 

5 

FRENCH  GERMAN  RESEARCH  INST 
ATTN  HANS-JURGEN  ERNST 

FRANCIS  JAMET 

PASCALE  LEHMANN 

K  HOOG 

H  LERR 

CEDEX  5  5  RUE  DU  GENERAL 

3 

CENTRE  D’ETUDES  DE  VAUJOURS 

ATTN  PLOTARD  JEAN-PAUL 

ERIC  BOTTET 

TAT  SIHN  VONG 

CASSAGNOU 

SAINT  LOUIS  68301 

FRANCE 

BOITE  POSTALE  NO  7 

COUNTRY  77181 

FRANCE 

1 

LABORATOIRE  DE  TECHNOLOGIE  DES 
SURFACES 

ECOLE  CENTRALE  DE  LYON 

ATTN  VINET  P 

6 

CENTRE  DE  RECHERCHES 

ET  D’ETUDES  D’ARCUEIL 

ATTN  D  BOUVART 

C  COTTENNOT 

BP  163 

69131  ECULLY  CEDEX 

FRANCE 

S  JONNEAUX 

H  ORSINI 

S  SERROR 

F  TARDIVAL 

16  BIS  AVENUE  PRIEUR  DE 

LA  COTE  D’OR 

1 

BATTELLE  INGENIEUTECHNIK  GMBH 
ATTN  W  FUCHE 

DUESSELDORFFER  STR  9 

ESCHBORN  D  65760 

GERMANY 

F941 14  ARCUEIL  CEDEX 

FRANCE 

1 

CONDAT 

ATTN  J  KIERMEIR 

MAXIMIL1ANSTR  28 

1 

DAT  ETBS  CETAM 

ATTN  CLAUDE  ALTMAYER 

ROUTE  DE  GUERRY  BOURGES 

8069  SCHEYERN  FERNHAG 

GERMANY 

18015 

FRANCE 

1 

DEUTSCHE  AEROSPACE  AG 

ATTN  MANFRED  HELD 

POSTFACH  13  40 

1 

ETBS  DSTI 

ATTN  P  BARNIER 

ROUTE  DE  GUERAY 

D  86523  SCHROBENHAUSEN 

GERMANY 

BOITE  POSTALE  712 

18015  BOURGES  CEDEX 

FRANCE 

1 

DIEHL  GBMH  AND  CO 

ATTN  M  SCHILDKNECHT 
FISCHBACHSTRASSE  16 

D  90552  ROTBENBACH  AD  PEGNITZ 

1 

FRENCH  GERMAN  RESEARCH  INST 

ATTN  CHANTERET  P-Y 

CEDEX  12  RUE  DE  IINDUSTRIE 

BP  301 

F68301  SAINT-LOUIS 

FRANCE 

GERMANY 
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4  ERNST  MACH  INSTITUT 
ATTN  VOLKER  HOHLER 
E  SCHMOLINSKE 
E  SCHNEIDER 
KTHOMA 
ECKERSTRASSE  4 
D-7800  FREIBURG  I  BR  791  4 
GERMANY 

1  EUROPEAN  SPACE  OPERATIONS  CENTRE 
ATTN  WALTER  FLURY 
ROBERT-BOSCH-STRASSE  5 

64293  DARMSTADT 
GERMANY 

3  FRAUNHOFER  INSTITUT  FUER 
KURZZEITDYNAMIK 
ERNST  MACH  INSTITUT 
ATTN  H  ROTHENHAEUSLER 
H  SENF 

E  STRASSBURGER 
HAUPTSTRASSE  18 
D79576  WEIL  AM  RHEIN 
GERMANY 

3  FRENCH  GERMAN  RESEARCH  INST 
ATTN  HARTMUTH  F  LEHR 
ROLF  HUNKLER 
ERICH  WOLLMANN 
POSTFACH  1260 
WEIL  AM  RHEIN  D-79574 
GERMANY 

2  IABG 

ATTN  M  BORRMANN 
H  G  DORSCH 
EINSTEINSTRASSE  20 
D  8012  OTTOBRUN  B  MUENCHEN 
GERMANY 

1  INGENIEURBURO  DEISENROTH 
AUF  DE  HARDT  33  35 
D5204  LOHMAR  1 
GERMANY 
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3  TU  CHEMNITZ-ZWICKAU 
ATTN  I  FABER 
L  KRUEGER 
LOTHAR  MEYER 

FAKULTAET  FUER  MASCHINENBAU  U. 

VERFAHRENSTECHNIK 
SCHEFFELSTRASSE  110 
09120  CHEMNITZ 
GERMANY 

1  TU  MUNCHEN 
ATTN  E  IGENBERGS 
ARCISSTRASSE  21 
8000  MUNCHEN  2 
GERMANY 

2  UNIVERSITAT  PADERBORN 
FACHBEREICH  PHYSIK 
ATTN  O  SCHULTE 

W  B  HOLZAPFEL 
D  33095  PADERBORN 
GERMANY 

1  BHABHA  ATOMIC  RESEARCH  CENTRE 
HIGH  PRESSURE  PHYSICS  DIVISION 
ATTN  N  SURESH 
TROMBAY  BOMBAY  400  085 
INDIA 

1  NATIONAL  GEOPHYSICAL  RESEARCH  INSTITUTE 

ATTN  G.  PARTHASARATHY 
HYDERABAD-500  007  (A.  P.) 

INDIA 

5  RAFAEL  BALLISTICS  CENTER 
ATTN  EREZ  DEKEL 
YEHUDA  PARTOM 
G  ROSENBERG 
Z  ROSENBERG 
Y  YESHURUN 
PO  BOX  2250 
HAIFA  31021 
ISRAEL 

1  TECHNION  INST  OF  TECH 
FACULTY  OF  MECH  ENGNG 
ATTN  SOL  BODNER 
TECHNION  CITY 
HAIFA  32000 
ISRAEL 
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1 

IHI  RESEARCH  INSTITUTE 

STRUCTURE  &  STRENGTH 

ATTN:  TADASHI  SHIBUE 

1-15,  TOYOSU  3 

KOTO,  TOKYO  135 

JAPAN 

3 

INSTITUTE  OF  CHEMICAL  PHYSICS 

RUSSIAN  ACADEMY  OF  SCIENCES 

ATTN  A  M  MOLODETS 

S  V  RAZORENOV 

A  V  UTKIN 

142432  CHERNOGOLOVKA 

MOSCOW  REGION 

1 

ESTEC  CS 

ATTN  DOUGLAS  CASWELL 

RUSSIAN  REPUBLIC 

BOX  200  NOORDWIJK 

2200  AG 

NETHERLANDS 

3 

INSTITUTE  OF  MECH  ENGINEERING  PROBLEMS 
ATTN  V  BULATOV 

D  INDEITSEV 

Y  MESCHERYAKOV 

2 

EUROPEAN  SPACE  AGENCY  ESTEC 

ATTN  LUCY  BERTHOUD 

MICHEL  LAMBERT 

POSTBUS  BOX  299  NOORDWIJK 

BOLSHOY,  61,  V.O. 

ST  PETERSBURG  199178 

RUSSIAN  REPUBLIC 

NL2200  AG 

NETHERLANDS 

1 

INSTITUTE  OF  MINEROLOGY  &  PETROGRAPHY 
ATTN  V  A  DREBUSHCHAK 

UNIVERSITETSKI  PROSPEKT,  3 

4 

PRINS  MAURITS  LABORATORY 

ATTN  H  J  REITSMA 

EDWARD  VAN  RIET 

630090  NOVOSIBIRSK 

RUSSIAN  REPUBLIC 

H  PASMAN 

R  YSSELSTEIN 

TNO  BOX  45 

RIJSWIJK  2280AA 

NETHERLANDS 

2 

IOFFE  PHYSICO  TECHNICAL  INSTITUTE 
DENSE  PLASMA  DYNAMICS 

LABORATORY 

ATTN  EDWARD  M  DROBYSHEVSKI 

A  KOZHUSHKO 

ST  PETERSBURG  194021 

1 

ROYAL  NETHERLANDS  ARMY 

ATTN  J  HOENEVELD 

RUSSIAN  REPUBLIC 

4 

V  D  BURCHLAAN  31 

PO  BOX  90822 

2509  LS  THE  HAGUE 

NETHERLANDS 

HIGH  ENERGY  DENSITY  RESEARCH  CTR 

1 

IPE  RAS 

ATTN  A  A  BOGOMAZ 

DVORTSOV AIA  NAB  18 

ST  PETERSBURG 

RUSSIAN  REPUBLIC 

ATTN  VLADIMIR  E  FORTOV 

GENADII I  KANEL 

V  A  SKVORTSOV 

O  YU  VOJOBIEV 

IZHORSKAJA  STR  13/19 

MOSCOW  127412 

2 

LAVRENTYEV  INST.  HYDRODYNAMICS 
ATTN  LEV  A  MERZHIEVSKY 

VICTOR  V  SILVESTROV 

630090  NOVOSIBIRSK 

RUSSIAN  REPUBLIC 

RUSSIAN  REPUBLIC 

1 

MOSCOW  INST  OF  PHYSICS  &  TECH 

ATTN  S  V  UTYUZHNIKOV 

1 

INSTITUTE  OF  CHEMICAL  PHYSICS 

ATTN  A  YU  DOLGOBORODOV 

KOSYGIN  ST  4  V  334 

MOSCOW 

RUSSIAN  REPUBLIC 

DEPT  OF  COMPUTATIONAL 

MATHEMATICS 

DOLGOPRUDNY  1471700 

RUSSIAN  REPUBLIC 
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1  RESEARCH  INSTITUTE  OF  MECHANICS 
NIZHNIY  NOVGOROD  STATE  UNIVERSITY 
ATTN  A  SADYRIN 
P.R.  GAYARINA  23  KORP  6 
NIZHNIY  NOVGOROD  603600 
RUSSIAN  REPUBLIC 

1  RUSSIAN  FEDERAL  NUCLEAR  CENTER 
ATTN  LEONID  F  GUDARENKO 
MIRA  AVE.,  37 
SAROV  607190 
RUSSIAN  REPUBLIC 

1  SAMARA  STATE  AEROSPACE  UNIV 
ATTN  L  G  LUKASHEV 
SAMARA 

RUSSIAN  REPUBLIC 

1  TOMSK  BRANCH  OF  THE  INSTITUTE 
FOR  STRUCTURAL  MACROKINETICS 
ATTN  V  GORELSKI 
8  LENIN  SQ  GSP  18 
TOMSK  634050 
RUSSIAN  REPUBLIC 

1  UNIVERSIDAD  DE  CANTABRIA 
FACULTAD  DE  CIENCIAS 
DEPARTMENTO  DE  FISICA  APLICADA 
ATTN  J  AMOROS 
AVDA  DE  LOS  CASTROS  S/N 
39005  SANTANDER 
SPAIN 

4  DEPARTMENTO  DE  QUIMICA  FISICA 
FACULTAD  DE  CIENCIAS  QUIMICAS 
UNIVERSIDAD  COMPLUTENSE  DE 
MADRID 

ATTN  V  G  BAONZA 
M  TARAVILLO 
M  CACERAS 
J  NUNEZ 
28040  MADRID 
SPAIN 
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1  CARLOS  III  UNIV  OF  MADRID 
ATTN  C  NAVARRO 
ESCUELA  POLITEENICA 
SUPERIOR 
a.  BUTARQUE  15 
28911  LEGANES  MADRID 
SPAIN 

1  UNIVERSIDAD  DE  OVIEDO 
FACULTAD  DE  QUIMICA 
DEPARTMENTO  DE  QUIMICA  FISICA  Y 
ANALITICA 
ATTN  E  FRANCISCO 
AVENIDA  JULIAN  CLAVERIA  S/N 
33006  -  OVIEDO 
SPAIN 

1  DYNAMEC  RESEARCH  AB 
ATTN  AKE  PERSSON 
P.O.  BOX  201 

S-151  23  SODERTALJE 
SWEDEN 

6  NATL  DEFENCE  RESEARCH  EST 
ATTN  LARS  HOLMBERG 
ULF  LINDEBERG 
LARS  GUNNAR  OLSSON 
L  HOLMBERG 
B  JANZON 
I  MELLGARD 
FOA  BOX  551 
TUMBA  S-14725 
SWEDEN 

2  SWEDISH  DEFENCE  RSCH  ESTAB 
DIVISION  OF  MATERIALS 
ATTN  S  J  SAVAGE 

J  ERIKSON 
STOCKHOLM  S- 17290 
SWEDEN 

1  K&W  THUN 
ATTN  W  LANZ 
ALLMENDSSTRASSE  86 
CH-3602  THUN 
SWITZERLAND 
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AWE 

ATTN  MICHAEL  GERMAN 
WAYNE  HARRISON 
FOULNESS  ESSEX  SS3  9XE 
UNITED  KINGDOM 

CENTURY  DYNAMICS  LTD 
ATTN  NIGEL  FRANCIS 
DYNAMICS  HOUSE 
HURST  RD 
HORSHAM 

WEST  SUSSEX  RH12  2DT 
UNITED  KINGDOM 

DERA 

ATTN  I  CULLIS 
FORT  HALSTEAD 
SEVENOAKS  KENT  TNI 4  7BJ 
UNITED  KINGDOM 


2  UNIVERSITY  OF  KENT 
PHYSICS  LABORATORY 
UNIT  FOR  SPACE  SCIENCES 
ATTN  PHILIPPE  GENTA 
PAUL  RATCLIFF 
CANTERBURY  KENT  CT2  7NR 
UNITED  KINGDOM 

7  INSTITUTE  FOR  PROBLEMS  IN 
MATERIALS  STRENGTH 
ATTN  S  FIRSTOV 
B  GALANOV 
O  GRIGORIEV 

V  KARTUZOV 

V  KOVTUN 

V  MILMAN 

V  TREFILOV 

3.  KRHYZHANO VSK Y  STR 
252142.  KIEV- 142 
UKRAINE 


DEFENCE  RESEARCH  AGENCY 
ATTN  W  A  J  CARSON 
I  CROUCH 
C  FREW 
T  HAWKINS 
B  JAMES 
B  SHRUBSALL 

CHOBHAM  LANE  CHERTSEY 
SURREY  KT16  OEE 
UNITED  KINGDOM 

UK  MINISTRY  OF  DEFENCE 
ATTN  GRAHAM  J  CAMBRAY 
CBDE  PORTON  DOWN  SALISBURY 
WITTSHIRE  SPR  OJQ 
UNITED  KINGDOM 


INSTITUTE  FOR  PROBLEMS  OF  STRENGTH 
ATTN  G  STEPANOV 
TIMIRYAZEVSKAYU  STR  2 
252014  KIEV 
UKRAINE 

INSTITUTE  OF  PHYSICS 

SILESIAN  TECHNICAL  UNIVERSITY 

ATTN  E  SOCZKIEWICZ 

44-100  GLIWICE 

UL.  KRZYWOUSTEGO  2 

POLAND 


K  TSEMBELIS 
SHOCK  PHYSICS  GROUP 
CAVENDISH  LABORATORY 
PHYSICS  &  CHEMISTRY  OF  SOLIDS 
UNIVERSITY  OF  CAMBRIDGE 
CAMBRIDGE  CB3  OHE 
UNITED  KINGDOM 
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